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Abstract 


An  investigation  was  made  to  determine  the  advantages 
of  modeling  the  wake  as  a  continuous  vortex  sheet  in  an 
unsteady  potential-flow  computer  model.  The  existing 
program  modeled  the  body  as  a  continuous  vortex  sheet,  and 
modeled  the  wake  as  discrete  vortex  cores  employing 
principles  of  vortex-lattice  methods.  A  method  is  developed 
to  approximate  the  wake  with  triangular  vortex  panels.  The 
method  accounts  for  redistributing  vortex  strength  as  the 
panels  expand  and  deform.  Comparisons  are  made  with  results 
of  the  original  program  for  rectangular  wings  at  varying 
angles  of  attack.  The  simulation  was  attempted  for  delta 
wings,  but  encountered  numerical  difficulties.  The  limited 
results  suggest  that  the  theory  and  method  are  valid,  and 
recommended  areas  are  given  for  further  pursuit. 
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MODELING  THE  WAKE  AS  A  CONTINUOUS  VORTEX  SHEET  IN 
A  POTENTIAL-FLOW  SOLUTION  USINC  VORTEX  PANELS 


I .  Introduction 


As  the  cost  of  producing  and  flight-testing  prototypes 
has  become  more  and  more  prohibitive,  the  increasing  power 
and  affordability  of  computers  make  them  an  invaluable 
alternative  in  the  engineering  design  process.  To  date,  a 
great  deal  of  experimental  data  is  available  and  is  used  to 
predict  the  aerodynamic  behavior  of  a  design,  but  a  robust 
mathematical  model,  particularly  in  the  realm  of  vortex- 
dominated  flows,  would  be  a  very  useful  tool.  The  present 
work  attempts  improvement  of  an  existing  method  developed  by 
Curtis  Mracek  (1)  in  1988.  He  developed  an  highly 
successful  unsteady,  potential-flow  aerodynamic  model  based 
on  vortex  panel  methods.  This  Fortran  program  models  the 
wake  development  by  way  of  discrete  vortex  cores  after  flow 
is  started  impulsively,  and  accurately  predicts  the  pressure 
distribution  in  vortex-dominated  flow  over  wings.  The 
present  research  pursues  Mracek' s  recommendation  to  model 
the  wake  as  a  continuous  vortex  sheet. 


Motivation  for  Research 


In  addition  to  the  financial  benefits  of  computer 

simulation  already  mentioned,  there  are  definite  time  and 

convenience  returns  to  be  exploited  in  the  design  process  if 

an  accurate  computer  model  exists.  Mracek's  model  was  a 

marked  advance  in  that  it  extended  previous  methods  limited 

to  steady  and  quasi-steady  flows  to  general  unsteady  flow, 

and  introduced  control  surface  motion  and  aerodynamic- 

dynamic  simulation.  Vortex-dominated  flows  are  of 

particular  interest  as  high  angle  of  attack  capabilities  are 

sought  in  modern  aircraft.  As  Mracek  states, 

Presently  there  are  several  good  methods  for 
predicting  steady  aerodynamic  loads  on  aircraft  at 
low  angles  of  attack  for  slow  speeds... The 
motivation  for  investigating  higher  angles  of 
attack  with  vortex  dominated  flow  is  that  many 
modern  aircraft  routinely  fly  in  this  regime,  and 
recently  efforts  have  been  made  to  harness  the 
vorticity-induced  forces  by  the  use  of  strakes, 
leading  edge  blowing,  and  vortex  channels.  Most 
numerical  methods  have  trouble  predicting  the 
forces  and  moments  for  these  conditions,  and  most 
cannot  be  used  for  unsteady  motion.  One  reason 
they  have  trouble  at  high  angles  of  attack  is  the 
flow  around  the  wing  is  heavily  influenced  by  the 
vortex  separation  along  the  leading  edge.  (1:2) 

The  existing  program  is  an  hybrid  vortex  method  in  that 

the  body  (wing)  is  discretized  into  triangular  panels  of 
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vorticity  that  varies  linearly,  and  the  wake  is  modeled  as  a 
mesh  of  discrete  vortex  cores  using  a  generalization  of  the 
procedure  employed  in  the  unsteady  vortex-lattice  method, 
will  be  shown  later,  the  vortex  panel  method  has  two  major 
advantages  over  using  discrete  vortex  cores  that  this 
research  attempts  to  exploit: 

(1)  The  disturbance  velocity  induced  by  a  vortex  panel 
approaches  infinity  at  a  slower  rate  as  the  distance 
from  the  panel  decreases  to  zero  since  it  is  a 
logrithmic  singularity.  Hence,  artificial  ''cutoff" 
distances  inside  which  induced  velocities  are 
arbitrarily  set  to  zero  can  be  reduced. 

(2)  The  disturbance  velocities  induced  by  a  vortex 
panel  are  more  smoothly  distributed  than  those  of 
discrete  cores  very  near  the  vorticity  itself.  Hence, 
nature  is  more  closely  approximated. 

Scope  of  Development 

Even  though  the  previously  existing  program 
incorporates  the  dynamics  of  moving  control  surfaces  and 
roll-,  pitch-,  and  yaw- rates,  this  development  considers 
only  impulsively  started  uniform  flow.  The  focus  is  on 
building  a  continuous  vortex  sheet  to  model  the  wake  similar 
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to  the  bound  vorticity  approximation  of  the  body.  This  is 
done  by  using  the  vorticities  already  calculated  by  the 
program  (and  convected  into  the  wake  as  discrete  cores) . 
Using  bound  vortex  panels  to  calculate  induced  velocities  is 
much  more  expensive  in  terms  of  computer  time  and  memory 
than  using  discrete  vortex  cores.  Therefore,  the  method 
developed  is  only  to  be  used  for  points  of  interest  near  the 
vorticity  itself  where  the  above  advantages  are  realized. 
This  necessitates  a  study  to  determine  "how  close  is  too 
close"  for  the  discrete  core  method. 
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II .  Aerodynamic  Mathematical  Model 


It  is  neither  the  intent  nor  within  the  scope  of  this 
report  to  provide  an  exhaustive  development  of  vortex-panel 
or  vortex-lattice  methods.  A  brief  summary  follows  of  the 
basic  principles  employed  by  the  existing  program. 

Governing  Equations 

The  program  is  a  potential-flow  solution  using  vortex 
distributions.  Using  the  development  of  Karamcheti  (2:254), 
the  following  field  equations  define  the  flow.  The 
continuity  equation  for  incompressible  flow  is 

div  9=0  (1) 

where  V  is  the  total  velocity  of  the  fluid.  The  total 
velocity  at  any  point  in  the  flow  is  the  sum  of  the 
freestream  velocity  and  the  disturbance  velocities  induced 
by  any  bodies  immersed  in  the  flow. 

V  =  V  +  v 

disturbance  freestream 

"  ^  +  Vfs  (2) 

The  vorticity  field,  n,  is  defined  as 

n  =  curl  V  =  curl  ( Vd  +  Vfs)  (3) 
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Since  the  freestream  is  irrotational ,  Equation  3  reduces  to 


n  =  curl  (4) 

A  third  vector  field,  A,  is  chosen  such  that 

V  =  curl  A  (5) 

The  continuity  equation  then  becomes 

div(curl  A)  =  0  (6) 

which  is  satisfied  for  any  A.  Now  the  problem  of  finding 
the  disturbance  velocity  reduces  to  finding  the  solution  to 

n  =  -v2A  (7) 

As  long  as  div  n  =  0,  Green’s  Theorem  gives 


where  fi  is  the  vcrticity  at  i,  r  is  the  point  of  interest, 
and  R  is  the  region  containing  the  vorticity.  Karamcheti 
(2:530-532)  establishes  that  if  the  vorticity,  n,  is 
contained  in  a  region  of  thickness,  e,  adjacent  to  the  body 
surface,  in  the  limiting  process  as  e  approaches  zero,  the 
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product  of  | n |  and  €  remains  constant  and  n  becomes  tangent 
to  the  surface.  Define 


7  (s)  = 


lim  | n | e 
n  |  -"oo,  e-»0 


(9) 


Then  Equation  (8)  can  be  rewritten  as  a  surface  integral 


A(?)  = 


4ir 


7  (5) 


r-s 


da 


(10) 


where  S  is  the  surface  of  the  body.  The  disturbance 
velocity  induced  by  a  continuous  vortex  sheet,  then,  is 


vd(?)  =  curl  A  =  -  curl 

4  n 


7(s) 


r-s 


da 


(11) 


Since  the  divergence  of  the  curl  of  a  vector  is 
identically  zero,  any  vector  field  that  represents  the 
surface  vorticity  will  satisfy  the  continuity  equation  as 
long  as 


div  7=0  (12) 

Two  other  conditions  that  must  be  satisfied  for  the 
flow  to  represent  reality.  They  are  the  no-penetration 
condition  at  the  body  surface 
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V  •  n  =  0  on  S  (13) 

where  n  is  the  normal  to  the  surface,  and  the  far-field 
condition 


Vd(«)  =  0  (14) 

When  a  thin  lifting  surface  is  immersed  in  the  flow, 
two  additional  conditions  must  be  met.  They  are  the 
conservation  of  circulation,  r,  in  the  wake  as  prescribed  by 
the  Theorem  of  Kelvin  and  Helmholtz 

DT 

-  =  0  (15) 

Dt 


and  the  Kutta  Condition 

ACp  =  0  (16) 

along  the  edge  where  the  wake  joins  the  wing. 

The  method  developed  by  Mracek  models  the  body  as  a 
collection  of  triangular  vortex  sheets,  the  corners  of  which 
lie  on  the  actual  surface  of  the  body.  The  vorticity 
distribution  across  the  triangle  is  a  linear  combination  of 
the  vorticities  at  each  corner,  or  "node.”  The  wake  is 
modeled  as  a  lattice  of  constant  strength  vortex  cores. 
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Chapter  2  of  Mracek  gives  the  full  development  of  the  matrix 
equations  which  are  solved  to  obtain  the  vorticity  at  each 
node  in  the  global  reference  frame.  The  system  of  equations 
is  overconstrained  and  is  solved  by  satisfying  the 
divergenceless  condition  (Equation  12)  exactly  and 
minimizing  the  error  in  the  no-penetration  condition 
(Equation  13)  via  a  least  square  method. 

The  Triangular  Element  of  Vorticitv  (l:14ff) 

The  present  method's  approach  in  modeling  the  wake  as 
a  continuous  vortex  sheet  is  to  parallel  the  methods  used  by 
Mracek  in  approximating  the  wing.  This  not  only  prevents 
"reinventing  the  wheel,"  but  also  allows  taking  advantage  of 
subroutines  already  in  place  and  validated  for  calculating 
induced  velocities.  For  this  reason,  a  summary  of  Mracek' s 
methods  follows. 

The  disturbance  velocity  given  by  Equation  11  is 
approximated  by  the  sum  of  the  velocities  induced  by  each  of 
the  triangular  elements.  Figure  1  depicts  one  such  element, 
all  of  which  have  a  local  x-axis  defined  along  its  longest 
edge  and  a  local  y-axis  defined  such  that  the  coordinates 
(b,c)  are  both  positive. 
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A 

i 

(b,0 

Y 

/ 

X 

X 

(0.0) 

(q,0) 

Figure  1  Triangular  Element  and  Local  Coordinate  System 
(1:14) 


The  surface  vorticity  vector,  7,  of  each  element  lies 
entirely  within  the  plane  of  the  element.  That  is,  in  the 
element's  local  frame, 


7  =  7X1  +  7yD  (17) 

The  vorticity  across  each  element  is  approximated  by  a 
linearly  varying  combination  of  the  vorticity  vectors  at 
each  node. 

7  =  ( 7xi  f  1  +  7x2^2  +7xjf3)I  +  (7^^  +  7y2f  2  +  7y3f3)5  (18) 

The  7xj  and  7yi  are  the  x-  and  y-components ,  respectively,  of 
7  at  corner  "i"  as  shown  in  Figure  2.  By  convention,  corner 
1  is  at  the  local  frame's  origin,  corner  2  is  at  (a,0),  and 
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corner  3  is  at  (b,c) .  The  fi  are  basis  functions  which 
equal  unity  at  their  respective  corners  and  zero  at  the 
other  corners. 

f,  =  a^  +  b,y  +  1 
f2  =  a2x  +  b2y 
f3  =  a3x  +  b3y 

where 


a,  =  -1/a 
a2  =  1/a 


a3  =  0 


b1  =  (b-a)/ac 
b2  =  -b/ac 
b3  -  1/c 


Figure  2  Triangular  Element  with  Vorticity 


(1:56) 
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The  divergence  constraint  (Equation  12)  can  now  be  written 


(21) 

Edge  Formulation  (1:39) 

For  there  to  be  non-zero  vorticity  perpendicular  to 
the  edge  of  a  continuous  vortex  sheet,  a  variable  strength 
vortex  filament  (core)  lies  along  the  edge  to  generate  or 
capture  the  vorticity.  The  circulation  (T)  around  this 
vortex  core  and  the  vorticity  strength,  7,  are  related  by 
(2:544) 

dr(x) 

7  ( X)  -  (22) 

dx 


37y  3  3 

div  7  =  —  +  —  =  2  a.7xj  +  E  bj7yi  =  0 

ax  ay  1=1  i=i 


where  T(x)  is  the  circulation  at  position  x  along  the  length 
of  the  core  and  7  is  the  strength  of  the  vortex  sheet.  The 
vorticity  is  perpendicular  the  edge  core  (i.e.,  7  =  7y  )  and 
varies  linearly  on  the  triangular  element  as  depicted  in 
Figure  3.  Along  the  edge. 


x  +  7y1 


(23) 
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Figure  3  Edge  Element  with  Associated  Edge  Core  (1:40) 


Integrating,  the  circulation  of  the  core  is 


7y2  "  7v1  ? 

r(x)  «  — 1 ( X2/ 2 )  +  7y1  x  +  G 
d  y1 

=  G,  (x2/  2)  +  G2  x  +  G3 


(24) 


G3  is  an  as  yet  unknown  constant  of  integration.  When 
modeling  a  body  (or  the  wake)  as  a  collection  of  triangular 
elements,  the  point  where  adjacent  edge  cores  meet  must 
exhibit  continuous  circulation  strength.  That  is,  for  two 
adjacent  elements,  n  and  n+1,  with  asssociated  edge  cores 
meeting  at  point  x0, 
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rn(Xo)  “  WX0>  (25) 

Because  of  the  relationship  between  circulation  and  vortex 
sheet  strength,  continuity  through  the  first  derivative  is 
also  enforced  so  that 


7y<n)(Xo)  7y<n+1)(Xo)  (26) 

These  constraints  actually  add  more  equations  than  unknown 
constants  of  integration  (G3's),  hence,  the  problem  remains 
overconstrained . 

Induced  Disturbance  Velocity 

The  model  of  a  body  immersed  in  the  flow  now  has  two 
components  —  the  triangular  elements  of  linearly  varying 
vorticity  and  the  quadratically  varying  vortex  filaments 
around  the  edges.  The  total  disturbance  velocity  induced  at 
any  point  in  the  flow  is  simply  the  sum  of  the  velocities 
induced  by  each  of  these  components. 

(27) 

vd  body  vd  elements  vd  edge  cores  '  ' 

The  velocity  induced  by  the  triangular  elements  is  given  by 
Equation  (11) . 
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The  velocity  induced  by  the  vortex  filaments  around 
the  edge  are  calculated  using  the  Biot-Savart  law.  For  each 
edge  filament,  a  local  axis  system  is  defined  with  the  x- 
axis  along  the  filament  in  the  direction  of  the  circulation. 
From  Equation  (24) ,  the  circulation  along  the  filament  is 
given  by 

F  (x)  =  G1  (x2/2)  +  G2  x  +  G3  (28) 

The  Biot-Savart  law  gives  the  velocity  induced  by  the 
filament  as 


- 


4  ir 


T(x) 


dl  x  (r-i) 
I r-i 1 3 


(29) 


where  dl  is  a  differential  length  along  the  filament 


dl  =  dx  l 


(30) 


when  expressed  in  the  local  coordinate  frame. 

The  Formation  of  the  Wake 

The  wake  is  a  vortex  sheet  emanating  from  the  edge  of 
the  lifting  surface.  It  must  satisfy  Equation  (12)  (the 
divergence  condition),  Equation  (14)  (far-field  condition), 
and  Equation  (15)  (conservation  of  circulation) .  The 
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existing  program  generates  the  wake  by  impulsively  starting 
the  flow  and  taking  snapshots  of  the  developing  wake  at 
discrete  time  intervals.  (1:65-69) 

It  must  be  understood  at  the  outset  that  the  intent  of 
the  present  study  is  not  to  change  the  method  of  generating 
and  convecting  the  wake.  Rather,  the  intent  here  is  to 
change  how  the  existing  wake  is  modeled  once  generated,  for 
the  purpose  of  calculating  velocities  induced  by  the 
vorticity  in  the  wake.  Any  changes  in  induced  velocities 
will,  of  course,  alter  the  development  of  the  wake,  but  the 
basic  principles  employed  in  developing  the  wake  remain  (1) 
the  generation  of  edge  cores  on  the  wing,  and  (2)  convecting 
them  downstream  at  the  local  particle  velocity. 

To  more  clearly  understand  the  focus  of  this  study, 
consider  a  wing  in  an  impulsively  started  flow.  At  the 
instant  the  flow  is  started,  vorticity  forms  on  the  surface, 
and  a  vortex  filament  (core)  with  circulation  develops 
around  the  edge  of  the  wing.  Not  enough  time  has  yet 
elapsed  to  move  this  circulation  off  the  wing  and  into  the 
flow.  After  an  incremental  time,  the  circulation  is  swept 
off  the  wing  and  downstream  in  the  flow,  and  a  vortex  sheet 
stretches  from  the  core  to  the  edge  from  whence  it 
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originated.  As  this  vorticity  moves,  its  location  relative 
to  the  wing,  hence,  the  flow  field  on  the  surface  of  the 
wing  is  changing.  At  the  end  of  the  first  time  increment, 
there  is  a  new  distribution  of  vorticity  on  the  surface,  so 
the  vortex  sheet  from  the  starting  core  to  the  surface  of 
the  wing  is  no  longer  compatible  with  the  edge.  To  resolve 
the  discontinuity,  a  new  edge  core  is  formed.  The  strength 
of  this  new  core  is  exactly  the  negative  of  the  strength  of 
the  starting  core  added  to  the  strength  of  the  core  needed 
to  generate  the  new  surface  vorticity  in  the  absence  of  any 
wake. 

At  the  start  of  the  second  time  step,  a  new  vortex 
distribution  is  formed  on  the  surface  because  of  the 
disturbance  induced  by  the  vorticity  in  the  wake.  Now  edge 
cores  are  again  formed  to  generate  and  capture  the  surface 
vorticity  as  in  the  first  time  step.  Additionally,  the 
negative  of  the  starting  vortex  is  collocated  at  the 
convecting  edge  to  capture  the  surface  vorticity  of  the  wake 
sheet.  Total  point  velocities  are  calculated  for  this 
configuration,  and  the  new  edge  cores  are  convected  into  the 
wake  along  with  the  starting  edge  core  and  the  vortex  sheet 
which  connects  them.  A  new  vortex  sheet  forms  connecting 


the  second  edge  core  to  the  edge  of  the  surface.  This 
procedure  is  followed  for  all  successive  time  steps. 

In  the  existing  program,  the  vorticity  in  the  wake  is 
not  modeled  as  a  continuous  vortex  sheet.  Rather,  the  wake 
sheet  is  discretized  into  vortex  cores  of  constant  strength. 
The  circulation  strength  chosen  for  this  discrete  core  is 
the  average  circulation  value  of  the  convected  variable 
strength  vortex  core.  A  four-sided  ring  of  this  circulation 
strength  and  composed  of  variable  strength  vortex  cores 
connects  the  starting  core  to  the  convecting  edge, 
satisfying  the  conservation  of  circulation  —  the  general 
unsteady  vortex  lattice  method. 

Figure  4  illustrates  the  formation  of  the  wake.  In 
Figure  4a,  the  flow  is  impulsively  started,  vorticity  forms 
on  the  surface,  and  the  associated  edge  cores  have 
circulation.  For  clarity,  only  the  trailing  edge  cores  are 
shown  to  be  convected,  and  they  are  shown  slightly  offset 
from  the  trailing  edge.  The  average  circulations  of  these 
edge  cores  are  then  calculated  and  convected  off  the 
trailing  edge.  A  four-sided  ring  of  circulation  is  formed 
connecting  the  starting  vortex  core  to  the  trailing  edge  and 
a  new  trailing  edge  core  is  formed  due  to  the  change  in 
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surface  vorticity  caused  by  the  presence  of  the  wake.  This 
condition  is  shown  in  Figure  4b.  The  four  corners  of  the 
ring  are  the  "nodes”  of  the  wake  whose  positions  are  updated 
each  time  step.  Figure  4c  shows  the  state  after  one  more 
time  step.  In  this  method,  only  the  average  circulations 
and  the  updated  positions  of  the  nodes  for  the  next  time 
step  are  recorded  to  fully  describe  the  wake. 

Given  the  circulation  of  the  variable  strength  core  in 
Equation  (28),  the  average  circulation,  ravg,  of  a  core  of 
length  d  is  found  by 

J0  Gy  (xz/2)  4-  G2  x  +  G3 
ravg  -  d 

=  G,  ( d2/ 6 )  +  G2  d  +  G3  (31) 

The  requirement  that  the  wake  be  single  valued  in  its 
pressure  distribution  (Equation  16)  is  satisfied  by 
convecting  at  the  local  particle  velocity.  For  each  of  the 
nodes  in  the  wake  and  along  the  trailing  edge,  the  local 
velocity  is  calculated  by 

V(t)  =  Vfreestream  +  Vd(b<xJy)  +  Vd(wake)  (32) 

and  their  new  positions  for  the  next  time  step  are 
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determined  by  the  first-order  difference  formula 

r (t+At)  =  r  (t)  +  V (t)  At  (33) 

where  At  is  the  time  increment  between  time  steps. 

The  thrust  of  this  work  lies  in  the  method  of 
calculating  Vd(wake)  in  Equation  (32)  above.  The  existing 
program,  as  stated  earlier,  uses  the  Biot-Savart  law  to 
determine  the  velocity  induced  by  each  of  the  constant 
strength  cores  of  Tavg  circulation  and  each  of  the  variable 
strength  filaments  which  form  the  four-sided  circulation 
rings.  Vd(Hake),  then,  is  simply  the  sum  of  these  induced 
velocities.  In  the  present  method,  with  the  wake  modeled  as 
a  continuous  vortex  sheet,  ^d(walce)  is  determined  in  the  same 
manner  as  Vd(body),  summing  the  velocities  induced  by 
triangular  elements  which  approximate  the  vortex  sheet,  and 
the  velocities  induced  by  edge  cores  surrounding  this  new 
"body."  The  following  example  will  not  only  illustrate  this 
method,  but  will  also  reveal  the  motivation  for  the  effort 
by  demonstrating  the  advantages  given  in  Chapter  1. 

An  Example  of  Equivalent  Vortex  Systems 

The  advantages  given  in  Chapter  1  of  using  vortex 
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panels  over  discrete  cores  in  determining  induced  velocities 
can  be  demonstrated  using  a  simple  contrived  example.  This 
example  was  used  by  Mracek  (1:72-74)  to  illustrate  the 
validity  of  using  discrete  vortex  cores  in  modeling  the 
wake. 

Consider  a  single  row  of  the  wake  composed  of  four 
square  elements  of  unit  dimension  as  depicted  in  Figure  5. 
For  simplicity,  the  entire  row,  modeled  as  a  continuous 
vortex  sheet,  is  planar  with  the  vortex  distribution  as 
shown.  For  each  segment,  a  local  coordinate  frame  is 
defined  with  the  positive  x-direction  in  the  direction  of 
the  circulation,  and  the  positive  y-direction  toward  the 
interior  of  the  element.  The  basic  equation  relating 
circulation  and  vortex  strength  is 


_d 

dx 


r(x)  =  -  7y 


(34) 


For  each  of  the  ten  discrete  cores,  the  vorticity 
distribution  is  written  in  its  local  coordinate  frame  and 
integrated  to  obtain  the  equation  for  r.(x).  For  segment 
number  1, 

7y1  =  -2x  -  1  (35a) 
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Figure  5  Continuous  Vortex  Sheet  Example 


(1:72) 


Integrating, 

r,(x)  =  2  ( x2/ 2 )  +  x  +  K,  (35b) 

where  K1  is  an  unknown  constant  of  integration.  Similarly, 


7y2  “X  “3/ 

=3 

r2(x)  =  1  ( X2/ 2 )  +  3X  +  k2 

(36) 

O, 

II 

X 

1 

r3(x)  =  - 1  ( X2/  2 )  +  4X  +  Kj 

(37) 

7y4  =  2X  -  3 

r4(x)  =  “2  ( X2/  2 )  +  3X  +  k4 

(38) 
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Note  that  for  r5(x),  since  all  vorticity  is  in  the  global  y 
direction,  in  the  local  coordinate  frame, 


fl 

o 

r5(x)  =  Kj 

(39) 

Continuing  around, 

7  y6  =  2X  +  1 

r6(x)  =  -  2  ( X2/  2 ) 

-  X 

+  K6 

(40) 

7y7  =  X  +  3 

r7(x)  =  -i(x2/2) 

-  3x  +  K7 

(41) 

7y8  =  -X  +  4 

r8(x)  =  i ( x2/ 2 )  - 

4x 

+ 

(42) 

7y9  -  “2X  +  3 

r9(x)  =  2 ( X2/ 2 )  - 

3x 

+  K, 

(43) 

■< 

o 

II 

O 

n 

o 

£ 

II 

X 

o 

(44) 

Now  the  circulation 

around  the  edge  must 

be 

continuous 

so 

that 

rio  (i)  =  r,(0) 

-  *10  =  K, 

(45) 

r,(i)  -  r2(0) 

=>  2  +  K,  =  K2 

(46) 

r2(i)  =  r3(0) 

=>  7/2  +  K2  =  K3 

(47) 

r3d)  =  r4(0) 

=>  7/2  +  K3  =  K4 

(48) 
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r4(i)  =  r5(0) 


2  +  k4  =  Ks 


(49) 


r5(i)  =  r6(o) 

=> 

*5  =  k6 

(50) 

r6(i)  =  r7(o) 

10 

+ 

j* 

O' 

11 

K7 

(51) 

r7(i)  =  r8(o) 

-7/2  +  K7 

=  Ka 

(52) 

r8(i)  =  r9(0) 

=> 

-7/2  +  Kg 

=  *9 

(53) 

r9d)  =  r10(O) 

I 

to 

+ 

II 

K10 

(54) 

This  yields  nine  independent  equations  in  the  ten  Kf 
unknowns.  To  get  the  tenth  equation,  impose  the  additional 
constraint  that 

K,  =  -K10  (55) 

Solving  this  system  of  linear  equations  for  Ki ,  the 
circulation  distributions  are  then 


r,(x)  =  X2  +  X  -  11/2 

(56) 

r2(x)  =  X2/2  +  3X  -  7/2 

(57) 

r3(x)  =  -X2/2  +  4X 

(58) 

r4(x)  =  -X2  +  3X  +  7/2 

(59) 
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r5(x)  =  11/2 


(60) 


r6(x)  =  -x2  -  x  +  11/2  (6i) 

r7(x)  =  -X2/2  -  3X  +  7/2  (63) 

rB(x)  =  X2/2  -  4X  (64) 

r9(x)  =  X2  -  3X  -  7/2  (65) 

ri0(x)  =  -11/2  (66) 

The  equivalent  discrete  vortex  core  model  can  now  be 
calculated  as  shown  in  Figure  6  where 

g,  =  =  f9  (67) 

g2  =  f2  =  f8  (68) 

g3  =  f3  =  f7  (69) 

g  4  =  =  r6  (70) 

with 

f.  =  /  r(  (x)  dx  (71) 

For  this  example,  the  length  of  each  segment  is  unity,  so 


26 


the  limits  of  integration  are  zero  to  one.  Solving  for  each 
f,-  (hence,  the  g,-'s) 

gi  =  “14/3,  g2  =  -11/6,  g3  =  11/6,  g4  =  14/3  (72) 

Figure  7  illustrates  the  equivalent  discrete  vortex  core 
system. 

The  discrete  core  and  continuous  sheet  representations 
are  compared  by  analysis  of  the  velocity  induced  in  a  planar 
region  normal  to  the  plane  of  the  vortex  sheet.  Figure  8 
shows  the  velocity  field  around  the  continuous  vortex  sheet 
as  calculated  by  the  existing  program  as  well  as  the 
velocity  field  for  the  equivalent  discrete  core  arrangement. 
Note  that  the  induced  velocity  fields  are  nearly  identical 
except  for  very  near  the  sheet  where  the  continuous  sheet 
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Figure  7  Equivalent  Constant  Strength  Vortex  Core  Arrangement 
(1:73) 


representation  gives  a  much  more  uniform  induced  velocity 
distribution.  This  demonstrates  two  of  the  statements  made 
in  Chapter  1  —  (1)  the  continuous  sheet  yields  induced 
velocities  more  closely  approximating  reality,  and  (2)  this 
advantage  is  only  realized  in  close  proximity  to  the  vortex 
sheet. 
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Figure  8  Induced  Velocity  Fields  of  Comparative  Example 
(1:74) 
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In  addition  to  illustrating  the  motivation  for  this 


study,  the  previous  example  also  revealed  that  the  existing 
program  already  has  much  of  the  capability  required  to 
handle  a  wake  modeled  as  a  continuous  vortex  sheet. 
Subroutines  are  already  in  place  which  can  calculate  the 
disturbance  velocities  induced  by  a  continuous  sheet,  given 
that  the  information  required  to  describe  that  continuous 
sheet  is  available.  After  implementation,  the  primary 
difference  between  the  program  of  this  effort  and  the  old 
version  is  in  the  method  of  determining  the  disturbance 
velocities  induced  by  the  wake.  The  new  version  sums  the 
influence  of  triangular  elements  which  approximate  the  wake 
and  the  edge  cores  with  circulation  that  surround  them,  just 
as  the  existing  program  does  for  the  wing. 

The  objective  of  this  effort,  then,  is  to  devise  a 
scheme  for  discretizing  the  wake  into  triangular  elements  of 
vorticity  that  approximate  the  wake  sheet  and  account  for 
the  variable  strength  cores  which  surround  it?  calculating 
and  recording  all  the  characteristics  required  by  the 
velocity  subroutines  used  for  the  wing.  The  algorithm  which 
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models  the  wing  serves  as  a  pattern  for  this  effort,  but 
modifications  must  be  made  since  the  wake  is  continually 
growing  and  distorting.  In  particular,  the  conservation  of 
circulation  is  maintained  in  the  existing  program  by  way  of 
the  constant  strength  r-average  values  which  convect  with 
the  four  sided  rings  of  circulation.  The  vortex  strength  of 
the  sheet  inside  that  ring  is  directly  related  to  that 
constant  strength  circulation.  The  first  obstacle,  then,  is 
to  somehow  distribute  the  vorticity  strength  over  an 
increasing  area  as  the  four  sided  ring  grows  with  time. 

Redistributing  the  Vorticitv  Over  an  Increasing  Area 

Figure  9  shows  the  convection  of  a  typical  edge  core  of 
length  d  off  the  trailing  edge  of  the  wing.  As  the 
endpoints  of  the  filament  are  convected,  their  updated 
positions  for  the  next  time  step  are  now  separated  by  a 
distance  of  d  +  Ad.  Associated  with  this  filament  are  the  G. 
coefficients  of  Equations  28  and  31  used  for  calculating 
ravg.  An  algorithm  must  be  devised  to  correct  these 
coefficients  for  the  change  in  filament  length  while 
maintaining  the  constant  Tavg. 
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That  is,  if  denoting  the  coefficients  corrected  for  growth 
as  primed, 

ravg  =  Gi  (d2/6)  +  G2(d/2)  +  G3 

=  G/ [  (d+Ad)2/6]  +  Gz'  [  (d+Ad)  /2  ]  +  G/  (74) 

The  first  attempt  was  to  simply  equate  the  terms  of 
like  powers  in  the  Tavg  equation.  This  yielded  the 
relationships 


=  G1 

[d/(d+Ad)]Z 

(75) 

=  G2 

[d/ (d+Ad) ] 

(76) 

"  G3 

(77) 

Upon  running  this  relationship,  it  was  found  that  there  were 


32 


marked  discontinuities  at  the  nodes  where  adjacent  filaments 
meet  which  violated  Equation  (25)  (circulation  strength  node 
compatibility)  and  Equation  (26)  (vorticity  strength  node 
compatibilty) .  The  method  made  no  consideration  for  node 
compatibility,  solving  for  each  filament  independently,  and, 
in  retrospect,  actually  forced  the  discontinuity  at  the 
nodes.  A  new  approach  was  taken  to  consider  the  entire  row 
of  the  wake  generated  in  a  time  step,  maintaining  the 
average  strength  and  node  compatibility  while  correcting  for 
growth  and  distortion. 

System  of  Linear  Equations  in  G-Primes.  To  solve  for 
the  coefficients  (G-primes)  corrected  for  wake  growth,  a 
system  of  linear  equations  in  G-primes  was  constructed  for  a 
row  of  the  wake  n  vortex  filaments  wide.  Figure  10  shows  a 
simple  row  of  only  n=2  elements  to  demonstrate  the 
conditions  enforced  in  this  system  of  equations.  The  first 
condition  is  that  ravg  remains  constant. 


=  r. 

1  avg  i  avg 


(78) 


The  second  and  third  conditions  enforce  node  compatibility 
in  both  T(x)  and  -r(x). 
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Figure  10  Two  Filament  Wake  Row 


r'j  (d+ad)  =  r'l+1(0) 

(79) 

7/i(d+ad)  =  7,j+1(0) 

(80) 

The  final  conditions  are  referred  to  as  "end  conditions." 
The  two  end  cores  are  constant  strength  (not  variable  in 
their  local  x)  such  that  for  the  left  end, 

rVo)  =  r,  (o)  =»  g'31  =  g31  (si) 

and  for  the  right  end  (the  n-th  filament) , 

r'n(d+Ad)  =  rn(d)  (82) 

For  each  convected  edge  core,  then,  there  are  three  g' ■ 


unknowns.  For  each,  the  available  equations  include  the 
ravg  equation  and,  except  for  the  last  core,  the  two  node 
compatibility  equations.  Together  with  the  two  end 
conditions,  a  system  of  3n  equations  has  been  constructed  in 
3n  unknowns. 

This  system  of  equations  is  represented  in  matrix  form 
in  the  program  so  that  a  standard  simultaneous  equation 
solver  could  be  employed.  By  loading  the  equations  into  the 
matrix  in  a  specific  order,  a  banded  matrix  results  with 
bandwidth  seven;  a  characteristic  which  can  be  exploited  by 
equation  solvers  to  decrease  computational  time.  The  first 
equation  loaded  into  the  matrix  is  the  left  end  condition 
(Equation  81) .  Then  stepping  through  each  filament  from 
left  to  right,  the  Tavg  equation  (Equation  78)  and  the  two 
node  compatibility  equations  (Equations  79  and  80)  are 
loaded.  For  the  last  (n-th)  filament,  only  the  Tavg  equation 
and  the  right  end  condition  (Equation  82)  are  loaded.  The 
resulting  coefficient  matrix  is  banded  to  ±  three  elements 
of  the  main  diagonal,  and  is  solved  by  a  subroutine  MSOLV 
(3)  which  takes  advantage  of  the  banded  characteristic.  For 
a  wake  row  n-filaments  wide,  the  matrix  equation  is 
constructed  as 
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(d+Ad) /(d+Ad) 


The  Fortran  algorithm  which  builds  and  solves  this 
matrix  equation  is  in  subroutine  CONVE  (Appendix  A) .  A  test 
case  was  run  with  a  rectangular  wing  with  unity  aspect  ratio 
at  5  degrees  angle  of  attack.  This  case  convects  twelve 
edge  cores  each  time  step,  four  from  each  side  and  four  from 
the  trailing  edge.  The  run  was  taken  to  60  time  steps  which 
ensured  a  steady  state  condition  in  the  wake.  The  G-prime 
coefficients  were  calculated  and  stored  using  Equation  (83) 
in  the  old  version  of  the  program,  but  they  were  not  used  in 
any  calculations  affecting  wake  convection  .  The  resulting 
coefficients  were  analyzed  by  plotting  the  core  strength 
(Equation  28)  and  the  vortex  sheet  strength  (Equation  34) 
for  the  trailing  edge  and  for  wake  rows  5,  10,  15,  and  20  of 
the  wake.  Figure  ll  demonstrates  that  the  core  strength 
remains  smooth  and  continuous  as  the  wake  spreads.  Figure 
12,  the  vortex  sheet  strengths,  is  even  more  revealing.  As 
the  wake  spreads,  not  only  is  node  compatibility  maintained, 
but  the  strength  of  the  sheet  is  being  redistributed  over 
the  increased  area  as  reflected  in  the  shallowing  of  the 
plot  as  the  wake  span  increases.  That  is,  the  area  under 
the  curve  is  remaining  constant.  The  odd  function  type 
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symmetry  about  the  midspan  indicates  that  spanwise  and 
directional  symmetry  is  also  maintained.  It  is  concluded 
that  the  method  to  redistribute  vortex  sheet  strength  as  the 
wake  expands  is  successful.  The  corrected  coefficients 
are  required  later  by  the  subroutines  which  calculate 
disturbance  velocities  induced  by  the  wake. 
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Figure  lla  Trailing  Edge  Vortex  Core  Strength 
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Discretizing  the  Wake 

Having  successfully  accounted  for  wake  growth  in 
redistributing  the  vortex  sheet  strength,  a  design  is 
developed  for  discretizing  the  wake  into  the  triangular 
elements  of  vorticity.  In  addition  to  expanding,  the  four 
sided  rings  of  circulation  are  not  planar  as  they  convect. 
In  fact,  they  may  become  quite  distorted,  especially  in  the 
outer  regions  of  the  span  where  the  wake  sheet  is  rolling 
up.  The  scheme  that  was  chosen  approximates  the  four  sided 
element  with  four  triangles  that  share  a  common  vertex  at 
the  "midpoint”  of  the  quadrilateral.  This  "midpoint"  is 
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simply  the  average  of  the  x-,  y- ,  and  z-coordinates  of  the 
four  corner  nodes.  It  is  envisioned  that  the  four  triangles 
will  provide  the  flexibility  required  to  adequately 
approximate  the  quadrilateral  even  in  cases  of  extreme 
twisting  as  may  be  encountered  in  the  rolling  up  of  the 
wake.  An  example  of  the  discretization  scheme  is  shown  in 
Figure  13,  depicting  two  wake  rows  four  edge  cores  wide, 
each  surrounded  by  its  own  edge  cores.  For  clarity,  this 
figure  shows  a  nonexistent  gap  between  the  two  rows. 
Hereafter,  the  upstream  side  of  the  row  is  referred  to  as 
the  "1-position"  which  correlates  with  subscripting  methods 
in  the  Fortran.  Similarly,  the  downstream  side  is  referred 
to  as  the  "2-position." 

In  order  to  properly  account  for  the  influence  of  the 
wake  when  calculating  induced  velocities,  each  row  of  the 
wake  must  be  treated  as  a  new  "body."  This  is  due  to  the 
fact  that  until  steady  state  is  reached,  the  strength  of  the 
surrounding  edge  cores  on  the  upstream  (1-position)  side  of 
a  given  row  is  not  the  same  as  the  strength  of  the  cores  on 
the  downstream  (2-position)  side  of  the  previous  wake  row 
(although  they  are  collocated) .  Therefore,  the  G-prime 
coefficient  matrix  equation  must  be  constructed  and  solved 
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Figure  13  Example  Wake  Row  Discretization 

separately  for  both  the  1-  and  2-position  filaments  for  each 
wake  row.  With  this  design  established,  all  that  remains  is 
to  provide  the  induced  velocity  calculating  subroutines  with 
the  characteristics  they  require  to  define  these  "additional 
bodies"  in  the  flow. 

Characteristics  of  the  Wake  Rows 

The  existing  subroutines  that  calculate  induced 
disturbance  velocities  for  a  continuous  vortex  sheet  require 
several  characteristic  values  which  define  its  geometry  and 
vorticity.  Once  again,  each  row  of  the  wake  is  treated  as 
an  individual  body  or  continuous  sheet.  Since  the  program 
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already  modeled  the  wake  as  discrete  cores,  much  of  the 
required  information  was  already  available  or  could  be 
readily  calculated  from  existing  values  to  make  the 
equivalent  continuous  vortex  system.  The  modeling  of  the 
continuous  system  looks  much  like  the  inverse  of  the  process 
in  the  example  at  the  end  of  Chapter  1,  which  transformed 
the  wake  into  the  equivalent  discrete  system.  Because  of 
the  nature  of  the  discretization  scheme,  some  assumptions 
and  simplifications  are  made. 

An  important  feature  of  this  work  is  that  the 
characteristics  of  the  wake  rows  are  appended  onto  the 
characteristic  arrays  defining  the  wing.  This  allows  using 
the  same  variable  names,  hence,  existing  subroutines  to 
calculate  induced  velocities.  Minimal  bookkeeping  effort 
(by  node  or  element  number  subscript)  is  required  to 
associate  characteristic  values  with  their  respective 
continuous  sheet.  The  wake  characteristics  and  any 
assumptions  or  simplifications  made  in  defining  them  follow. 
Except  where  otherwise  noted,  the  characteristic  values  are 
calculated  and  stored  in  subroutine  CONSHT  (Appendix  B) . 

Position  Coordinates.  The  updated  positions  of  each 
wake  node  were  already  calculated  in  subroutine  CONVE  in  the 
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existing  program.  These  are  determined  using  the  first- 
order  finite  difference  formula  given  in  Equation  (33) ,  and 
are  appended  each  time  step  to  the  arrays  that  define  the 
geometry  of  the  wing. 

Element  Orientation.  Several  arrays  define  the 
orientation  of  each  element.  A  direction  cosine  matrix 
relating  each  element’s  local  frame  to  the  global  frame  is 
computed  and  stored  using  an  algorithm  identical  to  that 
used  for  the  elements  that  approximate  the  wing.  For 
elements  which  have  one  side  on  an  outside  edge,  the 
direction  of  circulation  is  specified  for  the  edge  core 
along  that  side.  A  vector  normal  to  each  element  is 
calculated  by  cross  product  of  the  local  x-axis  with  the 
side  from  (a,0)  to  (b,c),  and  an  unit  normal  is  calculated 
for  each  node  by  averaging  the  normals  of  every  triangular 
element  that  shares  that  node.  Finally,  Mracek's  velocity 
subroutines  require  a  "c-factor"  (1:22-26),  calculated  using 
Equation  2.4-20  of  his  dissertation  (Ref.l).  The  c-factor 
is  a  measure  of  co-planarity  of  elements  sharing  a  given 
node,  and  is  calculated  for  each  vertex  of  each  triangle  as 


c  = 


r  2  ,  2  ,  2-,  1/2 

[wx  +w  +  wz  ] 


K2  +  wy2  +  wz2  +  (ut/nz)2]V2 


(84) 
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where  the  w's  are  the  x-,  y-,  and  z-components  of  the 
vorticity  vector  at  the  node  expressed  in  the  element's 
local  frame,  and  nz  is  the  z-component  of  the  average  normal 
at  the  node  expressed  in  the  element's  local  frame. 

Vorticitv.  Once  the  orientation  of  the  elements  is 
specified,  the  basis  functions  (Equations  19)  for  that 
element  are  computed.  The  vorticities  at  each  node  along 
the  1-  and  2-positions  are  calculated  using  Equation  (22) 
and  the  Gi  coefficients  which  have  been  corrected  for 
growth.  The  vorticities  are  calculated  in  the  edge  cores' 
local  frames  and  transformed  to  global  coordinates,  but  with 
one  modification.  Figure  14  shows  a  single  quadrilateral 
ring  in  the  wake.  The  vorticity  at  each  of  the  four  corners 
is  calculated  using  Equation  (22) ,  but  the  direction  is  not 
defined  perpendicular  to  the  edge  core  as  discussed  in 
Chapter  2.  Rather,  to  yield  more  accurate  induced 
velocities  directionally,  the  vorticity  is  rotated  into  the 
x-direction  of  the  left  and  right  triangular  elements  of  the 
quadrilateral,  more  closely  approximating  the  equivalent 
discrete  core  system.  Finally,  the  vorticity  vector  of  the 
newly  defined  "midpoint"  node  is  calculated  by  averaging  the 
global  components  of  the  four  corner  vorticities. 
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Figure  14  Rotation  of  Corner  Vorticity  Vectors 

Velocity  Induced  bv  the  Wake 

The  entire  wake  has  now  been  modeled  as  a  continuous 
vortex  sheet,  or  more  accurately,  as  several  rows  of 
continuous  vortex  sheets  (each  wake  row  is  a  continuous 
sheet) .  To  determine  the  particle  velocity  at  any  point  in 
the  flow,  the  contribution  of  the  wake  is  calculated  for 
inclusion  in  Equation  32.  The  existing  program  called  a 
subroutine  which  calculated  wake  induced  velocities  using 
the  discrete  core  vortex  system.  The  present  effort 
replaced  this  subroutine  with  one  nearly  identical  to  the 
routine  called  for  velocities  induced  by  the  body  (wing) . 
Subroutine  VELWK  is  included  in  Appendix  C. 


IV.  Results 


Modeling  the  wake  as  a  continuous  vortex  sheet  yielded 
mixed  results.  Several  computer  simulations  were  attempted 
for  unit  aspect  ratio  rectangular  and  delta  wings  at  5,  10, 
15,  and  20  degrees  angle  of  attack  for  comparison  with 
results  of  Mracek's  program.  The  simulations  encountered 
numerical  difficulties  for  the  rectangular  wing  at  20 
degrees  angle  of  attack  and  for  all  angles  of  attack  with 
the  delta  wing.  Successful  runs  were  made  with  the 
rectangular  wing,  however,  at  5,  10,  and  15  degrees.  These 
results  may  be  enough  for  analysis  of  the  theory  in  general 
by  comparing  with  results  of  the  program  that  models  the 
wake  as  discrete  cores.  It  appears  as  though  the  theory  may 
in  fact  be  valid  and  feasible  based  on  the  comparisons  as 
shown  in  Figures  15,  16,  and  17.  These  figures  show  cross 
sections  of  the  wake  at  several  constant  x-values  for  5,  10, 
and  15  degrees  angle  of  attack,  respectively.  Good 
correlation  is  demonstrated  for  these  limited  conditions, 
the  most  obvious  difference  being  the  amount  of  rolling  up 
in  the  wake.  The  suspected  causes  of  the  numerical 
difficulties  are  discussed  in  the  next  chapter. 
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Figure  15c  Wake  Cross  Sections  (A.O.A.  =  5  deg.,  x  =  -8) 
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Figure  16a  Wake  Cross  Sections  (A.O.A.  =  10  deg.,  x  =  -6) 
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Figure  16b  Wake  Cross  Sections  (A.O.A.  =  10  deg.,  x  =  -7) 
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Figure  16c  Wake  Cross  Sections  (A.O.A 


10  deg. ,  x  =  -8) 


Figure  16d  Wake  Cross  Sections  (A.O.A.  =  10  deg.,  x  =  -9) 
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Figure  17a  Wake  Cross  Sections  (A.O.A.  =  15  deg.,  x  =  -6) 


2.50  -i 
2.00 

1.50 
'■  .00 


o.co 


—  New  Progrom  a  =  15deg  x  = 

-  Old  Program  a=15deg  x  =  - 


-  / 
-  7 


'  j.  j  j  |  :  ; 

-4.00 


/  / 


Figure  17b  Wake  Cross  Sections  (A.O.A.  =  15  deg.,  x  =  -7) 
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Figure  17d  Wake  Cross  Sections  (A.O.A.  =  15  deg. ,  x  =  -9) 
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V.  Conclusions  and  Recommendations 


Conclusions 

An  algorithm  was  developed  to  model  the  wake  as  a 
continuous  vortex  sheet  in  a  potential-flow  solution  which 
had  modeled  the  wake  as  a  lattice  of  discrete  vortex  cores. 
Although  the  program  met  with  numerical  difficulties, 
results  obtained  from  limited  simulations  seem  to  indicate 
that  the  theory  and  method  are,  in  general,  valid  and  worthy 
of  further  pursuit.  Two  areas  are  suspect  as  being 
responsible  for  the  difficulties  encountered  at  high  angles 
of  attack  and  for  the  delta  wing. 

The  first  possible  shortcoming  is  that  the  present 
effort's  program  does  not  have  the  capability  to  "split" 
vortex  cores  in  the  wake  as  they  grow  beyond  a  pre- 
established  length.  The  original  program  has  this 
capability  which  smooths  the  approximation  of  the  wake, 
redistributes  the  vorticities,  and  allows  for  tighter 
rolling  up  of  the  wake  sheet.  The  exclusion  of  this 
capability  in  the  present  work  was  a  conscious  decision  for 
the  sake  of  simplifying  the  discretization  scheme  until  the 
theory  is  validated.  The  capability  could  be  re¬ 
established  by  revising  the  indexing  system  for  the  elements 


57 


and  nodes  in  the  wake.  This  implies  a  more  difficult 
bookkeeping  task,  but  not  an  adjustment  in  the  basic  theory.  . 

The  second,  and  probably  most  critical  area  is  the 
method  of  assigning  the  vorticity  at  the  newly  defined 
midpoints  of  the  guadrilateral  wake  elements.  By  simply 
assigning  the  average  of  the  four  corners,  the  divergence 
condition  of  Equation  (21)  is  not  satisfied  at  these  points. 
The  erroneous  vorticities  result  in  erroneous  induced 
velocities.  These  in  turn  distort  the  positions  and 
orientations  of  vorticities  for  the  next  time  step,  and  a 
snowball  effect  ensues.  This  difficulty  manifested  itself 
more  severely  in  high  angle  of  attack  and  delta  wing 
configurations  where  interaction  of  the  vortex  cores  is  more 
pronounced,  and  unrealistically  high  induced  velocities  were 
calculated. 

Both  of  these  shortcomings  may  be  contributing  in 
varying  degrees  to  the  error  in  the  simulation,  but  both  are 
reparable .  It  is  expected  that  their  correction  would 
eliminate  the  numerical  problems. 

Recommendations 

There  are  four  recommended  areas  for  further  study  in 
this  research.  The  first  two  are,  not  surprisingly, 


58 


concerned  with  the  two  problem  areas  discussed  above.  The 
most  critical  to  making  this  a  viable  tool  is  believed  to  be 
satisfying  the  divergence  condition  at  the  quadrilateral 
midpoints.  Reintroducing  the  wake-splitting  capability 
would  further  refine  the  simulation,  once  the  program  is 
operational . 

The  third  area,  mentioned  in  Chapter  1,  is  the  study  of 
"how  close  is  too  close"  for  the  discrete  core  method  to  be 
employed  in  calculating  induced  velocities.  This  effort, 
too,  is  secondary  to  resolving  the  numerical  difficulties, 
but  would  greatly  decrease  the  computer  time  required  for 
any  simulations.  This  would  require  developing  a  point-to- 
plane  distance  calculation  to  determine  whether  the  point  of 
interest  is  within  a  prescribed  cutoff  distance  of  the 
vortex  sheet. 

The  fourth  area  concerns  the  fact  that  this  program  is 
limited  to  symmetric  wings.  The  older  version  of  the 
program  had  this  limitation  due  to  the  method  of  assigning 
updated  position  coordinates  to  wake  nodes  for  the  next  time 
step.  A  similar  technique  for  mirroring  across  midspan  was 
used  in  this  work  for  assigning  vorticities  to  wake  nodes, 
further  entrenching  the  program  into  the  symmetric 
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requirement.  Minor  modifications  would  be  required  to 
remove  this  restriction,  making  it  a  more  powerful  design 
tool . 
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Appendix  A:  Subroutine  CONVE 


Subroutine  CONVE  convects  the  wake  by  calculating  the  total 
velocity  at  each  wake  position  and  determining  the  new  position  for  the 
next  time  step  using  the  first  order  finite  difference  formula.  The  Gi 
coefficients  corrected  for  expansion  and  deformation  of  the  wake  element 
are  calculated  for  both  the  1-  and  2 -positions  and  passed  along  with  the 
average  circulation  value  for  the  quadrilateral  wake  element.  To 
calculate  the  total  particle  velocity,  CONVE  calls  subroutine  VELLS 
(velocity  due  to  freestream)  and  VELBND  (velocity  due  to  the  bound 
vorticity  on  the  wing)  already  in  place  in  the  existing  program.  CONVE 
also  calls  VELWK  for  contribution  of  the  wake  (Appendix  C) . 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


************* 

*  CONVE  * 

************* 

THIS  SUBROUTINE  CONVECTS  THE  WAKE  AT  THE  LOCAL  VELOCITY. 
THE  Gi  COEFF.'s  OF  THE  G-AVERAGE  EQUATION  (CORRECTED  FOR 
ELEMENT  GROWTH)  ARE  CALCULATED  AND  CONVECTED  ALONG  WITH 
G-AVERAGES  AND  WAKE  POSITION  COORDINATES.  THIS  ROUTINE 
CALLS  SUBROUTINE  VELBND  TO  DETERMINE  VELOCITY  INDUCED  BY 
THE  BODY  AND  VELWK  TO  DETERMINE  VELOCITY  INDUCED  BY  THE 
WAKE  WHICH  HAS  BEEN  MODELED  AS  A  CONTINUOUS  VORTEX  SHEET. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  CONVE 
IMPLICIT  REAL*8 (A-H.O-Z) 

PARAMETER  (ME-1000 , MN-1000 , MEQ-200 , MUN-200 , MC-2 5 , MCI-100 ) 
PARAMETER  (MW- 2 5 , MWP-25 , MGP-240 , I DU- 7 , Ml-3 , M2-3 ) 

PARAMETER  (MWP1-101) 

DIMENSION  AGP(MGP ,7) , AGPT( 7 ,MGP) ,BGP(MGP) .MSH(MGP) ,GP(MGP) 
COMMON/VBOUND/  VBX,VBY,VBZ 
COMMON/VLSURF /  VLSX, VLSY, VLSZ 
COMMON/VEWAKE/  VWKX , VWKY, VWKZ 

C0MM0N/P01WAK/  X1WAK(MC ,MCI ,MW) . . . *AK(MC ,MCI ,MW) , Z1WAK(MC , MCI , MW) 
COMMON/ P02 WAK/  X2WAK ( MC , MC I , MW ) , Y2 WAK ( MC , MC I , MW ) , Z2WAK ( MC , MC I , MW ) 
COMMON/STRWAK/  GAVE(MC.MCI) ,CWAKE(MC ,MCI ,MW) , NRWAKE (MC , MW) 
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o  o  o 


COMMON/GNUMBE/ 

COMMON/WKCOEF/ 


COMMON /GBYEDG / 
COMMON/WKDIST/ 
COMMON/SPLITW/ 
COMMON/SPCOUN/ 
COMMON/SPCOU1/ 
COMMON/S P1WAK/ 
COMMON/S P2WAK/ 
COMMON / S  PGWAK / 
COMMON/SPWOR1/ 
COMMON/S PWOR2/ 
COMMON/EDDATA/ 
COMMON/EDGEDA/ 
COMMON/CONVEC/ 
COMMON/CONWOR/ 
COMMON/SP3WAK/ 
COMMON/VELINP/ 
COMMON/TIMESD/ 
COMMON/TIMES E/ 
COMMON/MODATA/ 
COMMON/ERRORS/ 
COMMON/CNTERS / 


G1(MC,MCI) ,G2(MC,MCI) ,G3(MC,MCI) 

G1WAKE(MC,MCI,MW) , G2WAKE(MC ,MCI ,MW) , 

G3WAKE(MC , MCI , MW) ,G1BAR(MC,MCI ,MW) , 
G2BAR(MC,MCI,MW) ,G3BAR(MC ,MCI ,MW) 

GIW(MC.MCI) ,G2W(MC,MCI) ,G3W(MC,MCI) 

DPDD(MC,MCI ,MW) ,D(MC,MCI ,MW) 

ISPLIT(MC,MCI,MW) , DIMINI , ISPER(MCI , MW) 

IATl(MW) , IAT2 (MW) , LAST1 (MW) ,LAST2(MW) 

NCOUNT(MW) .GCON(MCI) 

XIPER(MCI.MW) , Y1PER(MCI ,  MW) ,Z1PER(MCI,MW) 

X2PER(MCI ,MW) , Y2PER(MCI ,MW) , Z2PER(MCI ,MW) 
ISPL(MCI.MW) ,G(MCI ,MW) ,GPER(MCI ,MW) 

X1(MCI  fMW)  ,  Y1(MCI  ,MW)  ,Z1(MCI  ,MW) 

X2(MCI,MW) ,Y2(MCI,MW) ,Z2(MCI,MW) 

NEDGE ,NED(MCI , 2) , NCONV,NCO(MCI , 3) , NWKP , NWAKE1 (MWP) 
NOPEN , NCLOSE , NTOTCR , NEDG (MC , MCI , 2 ) , NCIRC (MC ) 
NCOVCR, NCONCR(MCI ) , NCON(MC ,MCI , 2) .NCONPO(MCI) 
XNP(MC ,MCI ,MWP1) , YNP(MC ,MCI ,MWP1) , ZNP(MC ,MCI ,MWP1) 
X2TEM(MC ,MCI ,MW) , Y2TEM(MC ,MCI ,MW) , Z2TEM(MC ,MCI ,MW) 
XP , YP , ZP , IDENTE , NODE1 , NODE2 

TIME , DTIME , NMOTIO , NCURTM , NTIME , NSTIME , NWPTS , MAXWK 
NSTMOT 

PMOMX , PMOMY , PMOMZ , SPAN , CHORD 
I ERROR , KERROR , MERROR , NERROR 
KCIR, NWELE 


OFF-O. 03D0*CHORD 


DOING  ONE  CIRCUIT  AT  A  TIME 

DO  900  KCIR-1, NCOVCR 

C 

C  INITIALIZING  THE  SPLITTING  OF  A  NEW  ROW 

C 

DO  20  I— 1 , NRWAKE(KCIR,NWPTS+1) -1 
ISPLIT(KCIR,I,NWPTS+l)-0 
20  CONTINUE 

C 

C  FINDING  THE  CONVECTED  WAKE  POSITIONS 

C 

DO  50  T  *1 .NWPTS+I 

C 

CCCCCCCC  THIS  IS  FOR  A  SYMMETRIC  WING 
C 

DO  30  J-l,NRWAKE(KCIR,I)/2+I 
XP-X1WAK(KCIR, J , I) 
YP-Y1WAK(KCIR,J, I) 
ZP-Z1WAK(KCIR,J,I) 

C 

C  FINDING  THE  TOTAL  VELOCITY 

C 


CALL  VELBND 
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n  o 


CALL  VELWK 
CALL  VELLS 
VAB SX-VBX+VWKX - VLSX 
VABSY— VBY+VWKY - VLSY 
VABS  Z-VBZ+VWKZ - VLS  Z 
C 

C  SAVING  THE  NEW  POSITIONS 

C 

XNP(KCIR,J,I)-XP+VABSX*DTIME 
YNP(KCIR,J,I)-YP+VABSY*DTIME 
ZNP ( KCIR , J , I ) — ZP+VABSZ*DTIME 
C 

C  SEEING  IF  THE  POINT  IS  TOO  CLOSE  TO  THE  WING 

C 

CALL  ISITOV(XNP(KCIR,J, I) ,YNP(KCIR,J , I) , ZNP(KCIR , J , I ) 
,OFF) 

30  CONTINUE 

DO  40  J-l ,NRWAKE(KCIR, I)/2 
K-NRWAKE ( KCIR ,  I )  -  J+l 
XNP(KCIR,K, I)-XNP(KCIR, J , I) 

YNP(KCIR,K, I)--YNP(KCIR, J , I) 

ZNP ( KCIR, K, I)-ZNP(KCIR, J , I) 

40  CONTINUE 

50  CONTINUE 

********************************************* 

IF(NWPTS.NE.O)THEN 

■A-****************************-**************** 

MOVING  THE  BAR  (1- POSITION)  CIRCULATION  BACK  TO  FRONT 

C 

ICO-O 

DO  320  II— 1 , NWPTS 
I-NWPTS-ICO 
IP1-I+1 

NSQ-NRWAKE ( KC I R , I ) - 1 
NEQ-3*NSQ 
c 

c  Create  "mshift"  array  for  "packing"  the  banded  matrix 
c  used  to  solve  for  G-primes. 

c 

DO  280  NK-l.NEQ 
MSH(NK)-0 

IF(NK  .GT.  4)  MSH(NK)-4-NK 
280  CONTINUE 
c 

c  Zeroing  out  the  matrix  equation  for  G-primes. 
c 

DO  285  ICL-l.NEQ 
BGP( ICL)-0 . 0D0 
DO  285  JCL-1,7 

AGP( ICL, JCL)=0 . 0D0 
285  CONTINUE 
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c 

c  Building  the  G-prime  matrix  equation  "packed".  The  equation 
c  is  banded  to  3  elements  either  side  of  the  main  diagonal.  The 
c  packing  offsets  the  rows  by  'mshift'  elements  to  form  a 
c  NEQ  by  7  matrix, 
c 

c  Load  row  one  of  G-prime  matrix  eq.--  the  left  end  condition, 
c 

AGP(1,3)— 1.0D0 
BGP(1)-G3BAR(KCIR, 1,1) 

DO  290  J-l.NSQ 
JP1-J+1 


Calculate  d  squared,  (d  +  delta  d)  squared  for  G-prime  Eq.'s. 


DSQ— ( (X1WAK(KCIR, JP1 , I) -X1WAK(KCIR, J , I) ) 
*(X1WAK(KCIR, JP1 , I) -X1WAK(KCIR, J , 1) ) )+ 
( (Y1WAK(KCIR, JP1, I) -Y1WAK(KCIR, J , I) ) 
*(Y1WAK(KCIR, JP1 , I) -Y1WAK(KCIR, J ,  I) ) )  + 
( (Z1WAK(KCIR, JP1 , I ) -Z1WAK(KCIR, J , I) ) 
*(Z1WAK(KCIR, JP1, I) -Z1WAK(KCIR, J , I) ) ) 
DPDDSQ— ( (XNP(KCIR, JP1 , I) -XNP(KCIR, J , I) ) 

*(XNP(KCIR, JP1 ,1) -XNP(KCIR, J ,  I) ) )  + 

( (YNP(KCIR, JP1 , I) -YNP(KCIR, J , I) ) 
*(YNP(KCIR, JP1 ,1) -YNP(KCIR, J , I) ) )+ 

( (ZNP(KCIR, JP1 , I) -ZNP(KCIR, J , I) ) 
*(ZNP(KCIR, JP1 , I) -ZNP(KCIR, J , I) ) ) 
DPDD(KCIR, J , IPl)-DSQRT(DPDDSQ) 

D(KCIR, J , I P 1 ) -DSQRT ( DSQ ) 


c  Indexing  to  build  banded  format 


IN1-3*J -1 
IN2-3*J -2 


c  Load  G- average  Eq.  into  matrix 
c 

AGP (INI , IN2+MSH(IN1) )-DPDDSQ/6 . 0D0 

AGP (INI , IN2+1+MSH(IN1) )-DPDD(KCIR, J , IPl)/2 .0D0 

AGP (INI , IN2+2+MSH(INl) )-l . 0D0 

BGP(INl)  — GWAKE(KCIR,  J  ,  I) 


IF  (J  . LT .  NSQ)  THEN 
c 

c  Load  Magnitude  Compatibility  Equation 
c 

IN1-IN1+1 

AGP ( INI , IN2+MSH(IN1) )-DPDDSQ/2 . 0D0 
AGP(IN1 , IN2+1+MSH( INI) )-DPDD(KCIR, J , IP1) 
AGP ( INI , IN2+2+MSH( INI) )-l . 0D0 
AGP( INI , IN2+5+MSH( INI) )-- 1 . 0D0 
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c 

c  Load  Slope  Compatibility  Eq. 
c 

IN1-IN1+1 

AGP(IN1,IN2+MSH(IN1))-DPDD(KCIR,J,IP1) 

AGP (INI , IN2+1+MSH(IN1) )-l .0D0 
AGP ( INI , IN2+4+MSH( INI) )— 1 . 0D0 
c 

c  If  last  element,  load  right  end  condition, 
c 

ELSE 

IN1-IN1+1 

AGP ( INI , IN2+MSH( INI ) ) -DPDDSQ/2 . 0D0 
AGP (INI , IN2+1+MSH( INI) )-DPDD(KCIR, J , IP1) 

AGP ( INI , IN2+2+MSH ( INI ) ) -1 . 0D0 

BGP(IN1)-G1BAR(KCIR, J , I)*DSQ/2 . 0D0+G2BAR(KCIR, J , I) 
*D(KCIR, J , IP1)+G3BAR(KCIR, J , I) 

END  IF 

290  CONTINUE 

c 

c  Transpose  AGP  for  solver 


DO  300  IT-1, NEQ 
DO  300  JT-1,7 

AGPT(JT,IT)-AGP(IT,JT) 

300  CONTINUE 

c  T 

c  Call  Banded  matrix  solver  for  [AGP]  [GP]-[BGP] 

c 

CALL  MSOLV( AGPT , GP , BGP , Ml , M2 , IDU , NEQ) 
c 

c  Assign  Position-1  coefficients  in  GBAR  arrays, 
c 

DO  310  J-l ,NRWAKE(KCIR, I) -1 
IND-3*J -  2 

G1BAR(KCIR, J , IPl)-GP(IND) 

G2BAR(KCIR, J , IP1)-GP(IND+1) 

G3BAR(KCIR, J , IP1)-GP( IND+2 ) 

310  CONTINUE 

ICO-ICO+1 
320  CONTINUE 
C 

C  MOVING  THE  CIRCULATION  BACK  TO  FRONT 

C 

ICO-O 

DO  120  II-l.NWPTS 
I-NWPTS -  ICO 
IP1-I+1 

NSQ-NRWAKE (KCIR, I) -1 
NEQ-3*NSQ 
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c 

c  Create  "mshift"  array  for  "packing"  the  banded  matrix 
c  used  to  solve  for  G-primes. 
c 

DO  80  NK-1 ,NEQ 
MSH(NK)-0 

IF(NK  .GT.  4)  MSH(NK)— 4-NK 
80  CONTINUE 
c 

c  Zeroing  out  the  matrix  equation  for  G-primes. 
c 

DO  85  ICL-l.NEQ 
BGP(ICL)-0.0D0 
DO  85  JCL-1,7 

AGP(ICL,JCL)-0.0D0 
85  CONTINUE 
c 

c  Building  the  G-prime  matrix  equation  "packed" .  The  equation 
c  is  banded  to  3  elements  either  side  of  the  main  diagonal.  The 
c  packing  offsets  the  rows  by  'mshift'  elements  to  form  a 
c  NEQ  by  7  matrix, 
c 

c  Load  row  one  of  G-prime  matrix  eq.--  the  left  end  condition, 
c 

AGP(1,3)-1.0D0 
BGP(l)— G3WAKE(KCIR, 1,1) 


DO  90  J-l.NSQ 
JP1-J+1 

c 

c  Calculate  d  squared,  (d  +  delta  d)  squared  for  G-prime  Eq.'s. 

c 

DSQ-( (X1WAK(KCIR, JP1 , IP1) -X1WAK(KCIR, J , IP1) ) 
*(X1WAK(KCIR, JP1, IP1) -X1WAK(KCIR, J , IP1) ) )  + 

( (Y1WAK(KCIR, JP1 , IP1) -Y1WAK(KCIR, J , IP1) ) 
*(Y1WAK(KCIR, JPl , IP1) -Y1WAK(KCIR, J , IP1) ) )+ 

( (Z1WAK(KCIR, JP1 , IP1) -Z1WAK(KCIR, J , IP1) ) 
*(Z1WAK(KCIR, JPl , IP1) -Z1WAK(KCIR, J , IP1) ) ) 
DPDDSQ— ( (XNP(KCIR, JPl , IP1) -XNP(KCIR, J , IP1) ) 

*(XNP(KCIR, JPl , IP1) -XNP(KCIR, J , IP1) ) )+ 

( (YNP(KCIR, JPl , IP1) -YNP(KCIR, J , IP1) ) 
*(YNP(KCIR , JPl , IP1) -YNP(KCIR, J , IP1) ) )+ 

( (JNP(KCIR, JPl , IP1) -ZNP(KCIR, J , IP1) ) 
*(ZNP(KCIR, JPl , IP1) -ZNP(KCIR, J , IP1) ) ) 
DPDD(KCIR, J , IPl)-DSQRT(DPDDSQ) 

D(KCIR, J , IPl)-DSQRT(DSQ) 
c 

c  Indexing  to  build  banded  format 
c 

IN1-3*J - 1 
IN2-3*J -  2 
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1 


Load  G- average  Eq.  into  matrix 

AGP ( INI , IN2+MSH( INI ) ) -DPDDSQ/6 . 0D0 

AGP( INI , IN2+1+MSH(IN1) )-DPDD(KCIR, J , IPl)/2 . 0D0 

AGP (INI , IN2+2+MSH ( INI ) )-l . 0D0 

BGP ( INI ) -GWAKE ( KCIR , J , I ) 


c 

c 

c 


c 

c 

c 


c 

c 

c 


IF  (J  .LT.  NSQ)  THEN 


Load  Magnitude  Compatibility  Equation 
IN1-IN1+1 

AGP ( INI , IN2+MSH ( INI ) ) -DPDDSQ/2 . 0D0 
AGP(IN1,IN2+1+MSH(IN1))-DPDD(KCIR,J,IP1) 
AGP (INI , IN2+2+MSH(INl) )-l . 0D0 
AGP (INI , IN2+5+MSH(INl) )--l .0D0 


Load  Slope  Compatibility  Eq. 

IN1-IN1+1 

AGP(IN1 , IN2+MSH(IN1) )— DPDD(KCIR, J , IP1) 
AGP (INI , IN2+1+MSH(IN1) )-l . 0D0 
AGP (INI , IN2+4+MSH ( INI ) )— 1 . 0D0 


If  last  element,  load  right  end  condition. 

ELSE 

IN1-IN1+1 

AGP (INI , IN2+MSH ( INI )) -DPDDSQ/2 . 0D0 

AGP ( INI , IN2+1+MSH ( INI ) ) -DPDD ( KCIR , J , I PI ) 

AGP ( INI , IN2+2+MSH ( INI ) ) -1 . 0D0 

BGP ( INI ) -G1WAKE ( KCIR , J , I ) *DSQ/2 . 0D0+G2WAKE ( KCIR , J , I ) 
*D(KCIR, J , IP1)+G3WAKE(KCIR, J , I) 

END  IF 


90  CONTINUE 

c 

c  Transpose  AGP  for  solver 
c 


100 


DO  100  IT-1, NEQ 
DO  100  JT-1,7 

AGPT (JT , IT ) -AGP ( IT , JT) 
CONTINUE 


c  T 

c  Call  Banded  matrix  solver  for  [AGP]  [ GP ] — [ BGP ] 

c 

CALL  MSOLV(AGPT,GP, BGP, Ml ,M2 , IDU,NEQ) 
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c  Assign  G-prime  answers  (Position-2  coefficients) 
c  to  appropriate  G#WAKE  arrays . 
c 

DO  110  J— 1 ,NRWAKE(KCIR, I) -1 
IND-3*J -  2 

G1WAKE(KCIR, J , IP1)— GP(IND) 

G2WAKE(KCIR,J,IP1)-GP(IND+1) 

G3WAKE(KCIR,J,IPl)-GP(IND+2) 

GWAKE(KCIR, J , IP1)-GWAKE(KCIR, J ,1) 

110  CONTINUE 

IC0-IC0+1 
120  CONTINUE 
END  IF 
C 

C  UPDATING  THE  POSITIONS  FIRST  THE  XI  POSITIONS 

C 

DO  70  I-l.NWPTS+l 
IP1-I+1 

DO  60  J-l , NRWAKE(KCIR , I) 

X1WAK(KCIR, J , IP1)-XNP(KCIR, J , I) 

Y1WAK(KCIR, J , IP1)-YNP(KCIR, J , I) 

Z1WAK(KCIR, J , IP1)-ZNP(KCIR , J , I) 

X2WAK(KCIR, J , I)-XNP(KCIR, J , I) 

Y2WAK(KCIR, J , I)-YNP(KCIR, J , I) 

Z2WAK(KCIR, J , I)-ZNP(KCIR, J , I) 

60  CONTINUE 

70  -CNTINUE 

NR WAKE (KCIR , NWPTS +2 ) -NRWAKE ( KC I R , NWPTS+1 ) 

*■*■***■*■*  >******************************************************** 
C  MOVING  THE  WAKE  STRENGTH  OFF  THE  WING  (G-prime  corrected 

C  for  expansion.) 

*****★•.  ********************************************************* 
NEQ— 3*NCONCR(KCIR) 
c 

c  Cre.te  "mshift"  array  for  "packing"  the  banded  matrix 
c  uset  to  solve  for  G-primes. 
c 

DO  130  NK-l.NEQ 
MSH(NK)-0 

IF(NK  .GT.  4)  MSri(NK)-4-NK 
130  CONTINUE 
c 

c  Zeroing  out  the  matrix  equation  for  G-primes. 
c 

DO  140  ICL-l.NEQ 
BGP( ICL)-0 . 0D0 
DO  140  JCL-1,7 

AGP(ICL,JCL)-0.0D0 
140  CONTINUE 
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c 

c 

c 

c 

c 

c 

c 

c 


Building  the  G-prime  matrix  equation  "packed".  The  equation 
is  banded  to  3  elements  either  side  of  the  main  diagonal.  The 
packing  offsets  the  rows  by  'mshift'  elements  to  form  a 
NEQ  by  7  matrix. 

Load  row  one  of  G-prime  matrix  eq.--  the  left  end  condition. 
AGP(1,3)— 1.0D0 

BGP(1)-G3(NC0N(KCIR, 1,1) ,NCON(KCIR, 1 , 2) ) 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


DO  150  J-l ,NC0NCR(KCIR) 
JP1-J+1 

L-NCON(KCIR, J , 1) 
M—NC0N(KCIR, J , 2) 


Calculate  d  squared,  (d  +  delta  d)  squared  for  G-prime  Eq.'s. 


DSQ-( (X1WAK(KCIR, JP1 ,1) -X1WAK(KCIR , J , 1) ) 
*(X1WAK(KCIR, JP1 , 1) -X1WAK(KCIR, J , 1) ) )  + 
( (Y1WAK(KCIR, JP1 , 1) - Y1WAK(KCIR, J , 1) ) 
*(Y1WAK(KCIR, JP1 , 1) -Y1WAK(KCIR, J , 1) ) )+ 
( (Z1WAK(KCIR, JP1, 1) -Z1WAK(KCIR, J , 1) ) 
*(Z1WAK(KCIR, JP1 , 1) -Z1WAK(KCIR, J , 1) ) ) 
DPDDSQ— ( (XNP(KCIR, JP1 , 1) -XNP(KCIR, J , 1) ) 

*(XNP(KCIR, JP1 , 1) -XNP(KCIR, J , 1) ) )+ 

( (YNP(KCIR, JP1 , 1) -YNP(KCIR, J , 1) ) 
*(YNP(KC1R, JP1 , 1) -YNP(KCIR, J , 1) ) )+ 

( (ZNP(KCIR, JP1 , 1) -ZNP(KCIR, J , 1) ) 
*(ZNP(KCIR, JP1 , 1) -ZNP(KCIR, J , 1) ) ) 
DPDD(KCIR, J , 1)— DSQRT(DPDDSQ) 

D(KCIR, J , l)-DSQRT(DSQ) 


Indexing  to  build  banded  format 

IN1-3*J - 1 
IN2-3*J-2 


Load  G-average  Eq.  into  matrix 

AGP (INI , IN2+MSH(IN1) )-DPDDSQ/6 . 0D0 

AGP (INI , IN2+1+MSH(IN1) )-DPDD(KCIR, J , l)/2 . 0D0 

AGP ( INI , IN2+2+MSH ( INI ) ) -1 . 0D0 

BGP( INl)-GAVE(L.M) 

IF  (J  .LT.  NCONCR(KCIR) )  THEN 


Load  Magnitude  Compatibility  Equation 
IN1-IN1+1 

AGP (INI , IN2+MSH(IN1) )-DPDDSQ/2 . 0D0 
AGP( INI , IN2+1+MSH( INI ) )-DPDD(KCIR , J , 1 ) 
AGP (INI , IN2+2+MSH( INI) )-l . 0D0 
AGP( INI , IN2+5+MSH( INI) )— 1 . 0D0 
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c  Load  Slope  Compatibility  Eq. 
c 

IN1-IN1+1 

AGP (INI , IN2+MSH( INI) )-DPDD(KCIR, J , 1) 

AGP (INI , IN2+1+MSH(IN1) )— 1 . 0D0 
AGP ( INI , IN2+4+MSH ( INI ) )  —  1 . 0D0 
c 

c  If  last  element,  load  right  end  condition, 
c 

ELSE 

IN1-IN1+1 

AGP ( INI , IN2+MSH ( INI ) ) -DPDDSQ/2 . 0D0 
AGP (INI , IN2+1+MSH(IN1) )-DPDD(KCIR, J , 1) 

AGP ( INI , IN2+2+MSH ( INI ) ) -1 . 0D0 

BGP(INl)-(Gl(L,M)*DSQ/2 . ODO)+(G2(L,M)*D(KCIR, J , 1) )+G3 (L,M) 
END  IF 

150  CONTINUE 
c 

c  Transpose  AGP  for  solver 
c 

DO  160  IT-1,  NEQ 
DO  160  JT-1,7 

AGPT(JT, IT)— AGP (IT, JT) 

160  CONTINUE 

c  T 

c  Call  Banded  matrix  solver  for  [AGP]  [GP]-[BGP] 

c 

CALL  MSOLV ( AGPT , GP , BGP , Ml , M2 , IDU , NEQ) 

c 

c  Assign  Position- 1  coefficients  in  GBAR  arrays  and  G-prime  answers 
c  (Position-2  coefficients)  to  appropriate  G#WAK£  arrays, 
c 

DO  170  J— 1 ,NRWAKE(KCIR, I) -1 
L-NC0N(KCIR, J , 1) 

M-NCON(KCIR, J , 2) 

G1BAR(KCIR,J,1)— G1(L,M) 

G2BAR(KCIR,J,1)— G2(L,M) 

G3BAR(KCIR, J , 1)--G3(L,M) 

IND-3*J -  2 

G1WAKE(KCIR,J, l)-GP(IND) 

G2WAKE(KCIR, J , 1)-GP(IND+1) 

G3WAKE(KCIR, J , l)-GP(IND+2) 

GWAKE(KCIR,J, l)-GAVE(L.M) 

170  CONTINUE 

900  CONTINUE 

RETURN 

END 
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Appendix  B :  Subroutine  CONSHT 


Subroutine  CONSHT  builds  the  wake  as  a  continuous  vortex  sheet  by 
appending  wake  information  onto  the  existing  arrays  that  define  the  body 
(wing)  as  a  bound  vortex  sheet.  These  values  are  required  by 
subroutines  VELWK  (Appendix  C) ,  and  VELE  and  VELVF  of  the  existing 
program  to  calculate  the  disturbance  velocity  induced  by  the  wake  when 
modeled  as  a  collection  of  triangular  elements  of  linearly  varying 
vorticity . 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

c  ****************  c 

C  *  CONSHT  *  C 

c  ****************  c 

c  c 

C  THIS  SUBROUTINE  BUILDS  THE  CONTINUOUS  WAKE  SHEET  BY  APPEND-  C 

C  ING  WAKE  INFORMATION  ONTO  ARRAYS  THAT  DEFINE  THE  BODY.  C 

C  C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  CONSHT 
IMPLICIT  REAL*8(A-H,0-Z) 

PARAMETER  (ME-1000 , MN-1000 , MEQ-200 , MUN-200 , MC-2  5 , MCI-100 ) 
PARAMETER  (MW-25 .MWP-25 ,MGP-240) 

COMMON/WKCOEF /  G1WAKE(MC ,MCI ,MW) ,G2WAKE(MC ,MCI ,MW) , 
G3WAKE(MC,MCI,MW) (G1BAR(MC,MCI,MW) , 
G2BAR(MC,MCI,MW) ,G3BAR(MC,MCI,MW) 

COMMON/GBYEDG /  G1W(MC,MCI) ,G2W(MC,MCI) ,G3W(MC,MCI) 

COMMON/ELDATA/  NELE , NODE , N(ME , 3 ) ,XNODE(MN) , YNODE(MN) , ZNODE(MN) 
COMMON/EDDATA/  NEDGE , NED (MCI , 2 ) , NCONV , NCO (MCI , 3 ) , NWKP , NWAKE1 (MWP ) 
COMMON/EDGEDA/  NOPEN, NCLOSE,NTOTCR,NEDG(MC, MCI , 2) ,NCIRC(MC) 
C0MM0N/P01WAK/  X1WAK(MC ,MCI ,MW) ,Y1WAK(MC ,MCI ,MW) , Z1WAK(MC ,MCI ,MW) 
C0MM0N/P02WAK/  X2WAK(MC , MCI , MW) , Y2WAK (MC , MCI , MW) , Z2WAK(MC , MCI , MW) 
COMMON/TIMESD/  TIME , DTIME , NMOTIO , NCURTM , NTIME , NSTIME , NWPTS , MAXWK 
COMMON/CONSTS/  FOURPI , PI , CUTOFF , CUTWK 

COMMON/CONVEC/  NCOVCR , NCONCR ( MCI ) , NCON (MC , MC I , 2 ) , NCONPO (MCI) 
COMMON/GNORMS/  ANORX(MN) , ANORY(MN) , ANORZ(MN) , CFACT (ME , 3 ) 
COMMON/STRWAK/  GAVE (MC, MCI) ,GWAKE(MC,MCI ,MW) , NRWAKE(MC.MW) 
COMMON/GOMEGB/  OMEX(MN) , OMEY ( MN ) ,OMEZ(MN) 

COMMON/GNUMBE/  Gl(MC.MCI) ,G2(MC,MCI) ,G3(MC,MCI) 

COMMON/CORNER/  AELE(ME) , BELE(ME) , CELE(ME) ,DIRCOS(ME,3 ,3) 
COMMON/EBASIS/  ABASE(ME, 6) 

COMMON/CONPOS/  XCP(ME) ,YCP(ME) , ZCP(ME) 

COMMON/WKDIST/  DPDD(MC , MCI , MW) , D(MC ,MCI , MW) 

COMMON/CNTERS/  KCIR.NWELE 


71 


*************************************************** 

C  Appending  the  x- ,  y- ,  and  z-location  of  wake 
C  nodes  onto  the  XNODE,  YNODE,  and  ZNODE  arrays. 
*************************************************** 

NNODE-NODE 

KCIR-1 

DO  40  I-l.NWPTS 

DO  10  J-l ,NRWAKE(KCIR, I) 

NNODE-NNODE+1 

XNODE (NNODE)— X1WAR (RCIR , J , I) 

YN0DE(NN0DE)-Y1WAR(RCIR, J , I) 

ZNODE (NNODE ) -Z1WAK(KCIR , J , I ) 

10  CONTINUE 
C 

C  Average  coords,  for  "midpoint"  of  4- sided  wake  element. 

C 

DO  20  J-l .NRWARE (RCIR, I) - 1 
JP1-J+1 
NNODE-NNODE+1 

XNODE ( NNODE ) - ( X1WAK ( KC IR , J , I ) +X2 WAR ( KC IR , J , I ) 

+X1WAR(RCIR, JP1, I)+X2WAR(RCIR, JP1 , I) )/4.0D0 
YN0DE(NN0DE)-(Y1WAK(KCIR, J , I)+Y2WAK(KCIR, J , I) 

+Y1WAR (RCIR , JP1 , I) +Y2WAR (KCIR , JP1 , I) )/4 . 0D0 
ZNODE(NNODE)— (Z1WAK(KCIR, J , I)+Z2WAK(KCIR, J , I) 

+Z1WAK(KCIR , JP1 , I )+Z2WAK(KCIR , JP1 , I) )/4 . 0D0 

20  CONTINUE 

DO  30  J—l .NRWARE (RCIR, I) 

NNODE-NNODE+1 

XN0DE(NN0DE)-X2WAK(KCIR, J ,1) 

YNODE (NNODE ) -Y2 WAR  <  KC IR , J , I ) 

ZNODE ( NNODE ) -Z2WAR ( RCIR , J , 1 ) 

30  CONTINUE 

40  CONTINUE 

************************************************************ 

C  Appending  the  wake  element  definitions  to  the  N(node,3) 

C  array  and  the  edge  data  for  each  wake  wake  row  to  NCIRC 

C  and  NEDG  arrays.  Also  create  Gltf,  G2W,  and  G3W  arrays 
C  of  wake  G-coeff.'s  subscripted  by  edge  core  #  for  VELBND 
************************************************************ 
NROW-NWPTS 
NWELE-NELE 
DO  60  I-l.NROW 
NBOD-NTOTCR+ I 

NC I RC ( NBOD ) - 2 *NRWARE ( RC I R , I ) 

NCOL-NRWARE ( RCIR , I ) - 1 

N1-2*NC0L+1 

N2-N1+1 

DO  50  J-l , NC0L 
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c 

1  *  -  * 

2 

• 

c 

c 

TRIANGLE  #1  OF  THE  WAKE  SQUARE 

\  / 

*  3 

NWELE-NWELE+1 

N ( NWELE , 1 ) -NOD E+ ( I - I ) * ( 3 *NCOL+ 2 ) + J 

N(NWELE , 2)-N(NWELE, 1)+1 

N (NWELE , 3 ) -NODE+ ( I - 1 ) * ( 3*NCOL+2 ) +NCOL+1+J 
NEDG (NBOD , J , 1 ) -N (NWELE , 1 ) 

NEDG ( NBOD , J , 2 ) -N ( NWELE , 2 ) 

GIW(NBOD, J)-G1BAR(KCIR, J , I) 
G2W(NBOD,J)-G2BAR(KCIR,J,I) 
G3W(NBOD,J)-G3BAR(KCIR,J,I) 

c 

* 

1 

c 

TRIANGLE  #2  OF  THE  WAKE  SQUARE 

3  *  | 

c 

NWELE-NWELE+1 

N (NWELE, 1)-N (NWELE -1,2) 

N (NWELE , 2 ) -NODE+ ( I - 1 ) * ( 3*NCOL+2 ) +2+NCOL+2+J 
N(NWELE, 3)-N(NWELE-l , 3) 

IF  (J  .EQ.  NCOL)  THEN 

* 

2 

NEDG ( NBOD , N2 , 1 ) -N ( NWELE , 2 ) 

NEDG ( NBOD , N2 , 2 ) -N ( NWELE , 1 ) 

G 1 W ( NBOD ,N2)-0.0D0 

G2W(NBOD , N2 ) — 0 . 0D0 

G3W(NBOD , N2)-G1WAKE(KCIR, J , I)*DPDD(KCIR, J , I)*DPDD(KCIR, J , I)  1 

/2 . 0D0+G2WAKE(KCIR , J , I) *DPDD (KCIR , J , I) 

• 

c 

+G3WAKE(KCIR, J , I) 

END  IF 

*  3 

c 

TRIANGLE  #3  OF  THE  WAKE  SQUARE 

/  \ 

c 

NWELE-NWELE+1 

N (NWELE , 1 ) -N (NWELE - 1 , 2 ) 

N ( NWELE , 2 ) -NODE+ ( I - 1 ) * ( 3*NCOL+2 ) +2*NCOL+l+J 

N (NWELE, 3)-N(NWELE- 1,3) 

NEDG (NBOD , J+NCOL , 1 ) -N( NWELE , 2 ) 

NEDG (NBOD, J+NCOL, 2) -N (NWELE, 1) 

G1W(NB0D , J+NCOL) -G1WAKE(KCIR , J , I) 

G2W(NBOD , J+NCOL) -G2WAKE( KCIR , J , I) 

G  3 W ( NBOD , J+NCOL) -G  3 WAKE ( KC I R , J , I ) 

2  *  -  * 

1 

c 

2  * 

c 

TRIANGLE  #4  OF  THE  WAKE  SQUARE 

1  * 

3 

c 

NWELE-NWELE+1 

N (NWELE , 1)-N(NWELE-1 , 2) 

N(NWELE, 2)-N(NWELE-3 , 1) 

N (NWELE, 3)-N(NWELE-l , 3) 

IF  (J  .EQ.  1)  THEN 

1  * 

NEDG ( NBOD , N1 , 1 ) -N ( NWELE , 2 ) 

NEDG ( NBOD, Nl, 2 )-N( NWELE, 1) 

• 
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GlW (NBOD , N1 ) —0 . 0D0 
G  2W ( NBOD , N 1 ) — 0 . 0D0 
G3W(NB0D,N1)— G3WAKE(KCIR, J , I) 

END  IF 

50  CONTINUE 
60  CONTINUE 

FINDING  THE  LONGEST  SIDE  FOR  EACH  ELEMENT 
DO  70  I-NELE+1 , NWELE 

PUTTING  IN  THE  APPROPRIATE  NODE  POSITIONS  INTO  INTERMEDIATE 
VARIABLES 

XI— XNODE(N(I , 1) ) 

Yl— YNODE(N(I , 1) ) 

Zl— ZNODE(N(I , 1) ) 

X2— XNODE(N(1 , 2) ) 

Y2-YNODE(N(1 , 2) ) 

Z2— ZNODE(N(I , 2) ) 

X3-XNODE(N(I , 3) ) 

Y3-YNODE(N(I , 3) ) 

Z3— ZNODE(N(I , 3) ) 

ROTATING  THROUGH  THE  POSITIONS  TO  FIND  THE  LONGEST  SIDE 

R1-(X2-X1)*(X2-X1) 

+(Y2-Y1)*(Y2-Y1) 

+(Z2-Z1)*(Z2-Z1) 

R2-(X3-X2)*(X3-X2) 

+(Y3-Y2)*(Y3-Y2) 

+(Z3-Z2)*(Z3-Z2) 

R3-(X1-X3)*(X1-X3) 

+(Y1-Y3)*(Y1-Y3) 

+(Z1-Z3)*(Z1-Z3) 

IF(R1 . LT . R2 . OR.R1 . LT . R3)  THEN 

IF(R2 . GT . R1 . AND . R2 . GT . R3 )  THEN 

FOR  THE  CASE  WHERE  THE  SECOND  VECTOR  IS  LONGEST 

NDUM-N (1,1) 

N(I,I)-N(I,2) 

N(I,2)-N(I,3) 

N(I , 3)— NDUM 
ELSE 

IF(R3 . GT . R1 . AND . R3 . GT . R2)THEN 

FOR  THE  CASE  WHERE  THE  THIRD  VECTOR  IS  LONGEST 
C 

NDUM-N (1,1) 

N ( I , 1 ) — N (1,3) 

N(I , 3)-N(I ,2) 

N ( I , 2 ) -NDUM 
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1 


END  IF 
END  IF 
END  IF 
70  CONTINUE 

***************************************** 

C  Appending  direction  cosine  matrices 
C  for  wake  elements  to  DIRCOS  array. 
***************************************** 

DO  80  I-NELE+1 , NWELE 
C 

PUTTING  IN  THE  APPROPRIATE  NODE  POSITIONS  INTO  INTERMEDIATE 
VARIABLES 
C 

Xl-XNODE(N(I,l)) 

Yl«YNODE(N(I , 1) ) 

Zl— ZNODE(N(I , 1) ) 

X2-XN0DE(N(I , 2) ) 

Y2-YNODE(N(I , 2) ) 

Z2— ZNODE(N(I , 2) ) 

X3— XNODE (N (1,3)) 

Y3— YNODE(N(I , 3) ) 

Z3-ZNODE (N ( 1,3)) 

C 

C  DEFINING  THE  INFLUENCE  MATRIX 

C 

D1-DSQRT((X2-X1)*(X2-X1)+(Y2-Y1)*(Y2-Y1)+(Z2-Z1)*(Z2-Z1)) 

AN1-(Y2-Y1)*(Z3-Z2)-(Z2-Z1)*(Y3-Y2) 

AN2-(Z2-Z1)*(X3-X2) - (X2-X1)*(Z3-Z2) 

AN3-(X2-X1)*(Y3-Y2)  -  (Y2-Y1)*(X3-X2) 

D3-DSQRT (AN1*AN1+AN2*AN2+AN3*AN3 ) 

DIRCOS (I , 1 , 1)-(X2-X1)/D1 
DIRC0S(I,1,2)-(Y2-Y1)/D1 
DIRCOS (I , 1, 3)— (Z2-Z1)/D1 
DIRCOS (I , 3 , 1)-AN1/D3 
DIRCOS (I , 3 , 2)-AN2/D3 
DIRCOS (I,3,3)-AN3/D3 

DIRCOS (1,2, 1) -DIRCOS (1,3,2) *DIRC0S (1,1,3) 

- DIRCOS (1,3,3) *DIRCOS (1,1,2) 

DIRCOS (I , 2 , 2) -DIRCOS (I , 3 , 3)*DIRCOS(I ,1,1) 

-DIRCOS (I , 3 , l)*DIRCOS(I ,1,3) 

DIRCOS (1,2,3) -DIRCOS (1,3,1) *DIRC0S (1,1,2) 

- DIRCOS ( I , 3 , 2 )*DIRCOS ( I , 1 , 1 ) 


C 

C  SETTING  THE  CONTROL  POINT  LOCATION  USING  THE  MIDPOINT  OF  THE 

C  ELEMENT 

C 

XCP(I)-(X1+X2+X3)/3.0D0 
YCP(I)-(Yl+Y2+Y3)/3 . 0D0 
ZCP(I)-(Z1+Z2+Z3)/3,ODO 
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c 

C  FINDING  THE  VERTICIES  OF  THE  TRIANGLE  NAMELY  A,B,C 

C 

AELE(I)— (X2-Xl)*DIRCOS(I ,1,1) 

+(Y2-Yl)*DIRCOS (1,1,2) 

+ ( Z2 - Z1 ) *DIRCOS (1,1,3) 
BELE(I)-(X3-Xl)*DIRCOS (I ,1,1) 

+(Y3-Yl)*DIRCOS(I ,1,2) 

+(Z3-Zl)*DIRCOS(I ,1,3) 
CELE(I)-(X3-Xl)*DIRCOS (I ,2,1) 

+(Y3-Yl)*DIRCOS (I , 2,2) 

+(Z3-Zl)*DIRCOS (I ,2,3) 

C****************************************** 

C  Appending  the  wake  elements'  basis 
C  function  coefficients  to  ABASE 

C** -A***********  jt*************************** 

ABASE (I , 1)— 1.0D0/AELE(I) 

ABASE ( I , 2 ) -1 . ODO/AELE ( I ) 

ABASE ( I , 3 ) — 0 . 0D0 

ABASE (I ,4)— (BELE(I) -AELE(I) ) /(AELE(I)*CELE(I) ) 
ABASE ( I , 5 )— -BELE( I) /(AELE( I )*CELE( I ) ) 

ABASE(I , 6)— 1 . 0D0/CELE(I) 

80  CONTINUE 

C****************************************** 

C  Appending  the  average  normals  for  wake 
C  nodes  to  ANORX.Y,  and  Z  arrays. 

C**  ******************************************  ■** 

DO  100  I— NODE+1 , NNODE 
SUMX-0 . 0D0 
SUMY-0 . 0D0 
SUMZ-0 . 0D0 
KOUNT-O 

DO  90  J-NELE+1 , NWELE 

C 

C  SEEING  IF  THE  ELEMENT  HAS  THE  NODE 

C 

IF(N(J,1) .EQ.I.OR. 

N(J,2) .EQ.I.OR. 

N(J , 3) . EQ. I) THEN 
C 

C  ADDING  IN  THIS  ELEMENTS  NORMAL 

C 

SUMX-SUMX+DIRCOS ( J . 3 , 1 ) 
SUMY-SUMY+DIRCOS ( J ,3,2) 
SUMZ-SUMZ+DIRCOS (J , 3 , 3 ) 

KOUNT -KOUNT  + 1 
END  IF 

90  CONTINUE 

ANORX ( I ) -SUMX/FLOAT ( KOUNT ) 

ANORY ( I ) -SUMY/FLOAT ( KOUNT ) 

ANORZ( I ) -SUMZ/FLOAT ( KOUNT ) 

100  CONTINUE 


76 


c 

C  MAKING  UNIT  NORMALS  AT  THE  NODES 

C 

DO  110  I-NODE+1 , NNODE 

AMAG-DSQRT ( ANORX ( I ) *ANORX ( I ) 

+ANORY ( I ) * ANORY ( I ) 

+ANORZ ( I ) *AN0R2 ( I ) ) 

ANORX ( I ) -ANORX ( I ) /AMAG 
ANORY ( I ) -ANORY ( I ) /AMAG 
ANORZ ( I ) -ANORZ ( I ) /AMAG 
110  CONTINUE 

C****************************************** 

C  Appending  wake  node  vorticities  onto 
C  the  OMEX.Y,  and  Z  arrays. 

C****************************************** 

c 

c  Zero  out  the  arrays 
c 

DO  130  I-NODE+1 , NNODES 
OMEX(I)— 0.0D0 
OMEY  ( I )  — 0 . 0D0 
OMEZ(I)-0 . 0D0 
130  CONTINUE 

DO  150  I— 1 , NWPTS 

NCOL-NRWAKE ( KC IR , I ) - 1 
DO  140  J-l.NCOL/2 

C 

C  The  "1- Positions" 

C  **************** 

C  Determine  the  element  #  (triangle  1  of  wake  square) 

C 

LEL-NELE+4*NC0L* ( I - 1 ) +4* ( J - 1 ) +1 
c 

c  Determine  the  node  number  for  local  (a,0) 
c 

IND-NODE+ ( I - 1 ) * ( 3*NCOL+2 ) +J 
c 

c  The  vorticity  strength  at  x-a.  The  G-coefficint 
c  arrays  define  the  strength  in  the  negative  local 
c  y-direction  for  triangle  1  of  the  squares  with  0 
c  at  the  (a,0)  position.  Convert  this  strength  to 
c  global  coords  in  the  positive  x-direction  of 
c  triangle  4  of  the  square, 
c 

GKF"  •G2BAR(KCIR,J,I) 

OMEA ( IND ) -GKEY+DIRCOS ( LEL+  3 ,1,1) 

0MEY( IND) -GKEY*DIRC0S ( LEL+3 ,1,2) 

OMEZ ( IND ) -GKEY+DIRCOS ( LEL+3 ,1,3) 


--*(0,0) 

/ 

c) 


77 


If  at  mid- span,  convert  the  strength  at  local  (0,0)  to  global 
coords  in  the  negative  x-direction  of  triangle  2  of  the  wake 
square.  The  strength  is  equal  to  -G2bar  of  the  next  wake  square. 


IF(J  .EQ.  NCOL/2)  THEN 
GKEY—  G2BAR(KCIR ,  J+l ,  I) 

OMEX ( IND+1 )-GKEY*DIRCOS (LEL+1 ,1,1) 
OMEY ( IND+1 ) -GKEY*DIRC0S (LEL+1 ,1,2) 
0ME2 ( IND+1 ) -GKEY*DIRC0S (LEL+1 ,1,3) 
IF(M0D(NC0L, 2)  .EQ.  0)  THEN 
0MEX( IND+1 )-0.0D0 
OMEY (IND+1) -0.0 DO 
0MEZ( IND+1 )-0.0D0 
END  IF 
END  IF 


The  "2-Positions" 
**************** 

Determine  the  element  # 
(triangle  3  of  wake  square) 


(0,0) 


(b,c) 

* 

/  \ 


*  (a,0) 


MEL-NELE+4+NC0L* ( I - 1 ) +4* ( J - 1 ) +  3 


Determine  the  node  #  for  local  x-0 

INN-N0DE+ ( I - 1 ) * ( 3*NCOL+2 ) +2*NC0L+1+J 


Convert  the  strength  at  local  (0,0)  to  global  coords  in  the 
negative  x-direction  of  triangle  4  of  the  wake  square. 

0MEX ( INN ) - - G2WAKE ( KC IR , J , I) +DIRC0S (MEL+1 ,1,1) 
0MEY(INN)--G2WAKE(KCIR, J , I)*DIRC0S(MEL+1 ,1,2) 
0MEZ(INN)--G2WAKE(KCIR, J , I)*DIRC0S(MEL+1 ,1,3) 


If  at  mid-span,  convert  the  strength  at  local  (a,0)  to  global 
coords  in  the  positive  x-direction  of  triangle  2  of  the  wake  square. 
The  strength  is  equal  to  G2- prime  of  the  next  wake  square. 


IF(J  .EQ.  NC0L/2)  THEN 

0MEX(INN+1)— -G2WAKE(KCIR, J+l , I)*DIRC0S (MEL- 1,1,1) 
0MEY( INN+1 )--G2WAKE(KCIR , J+l , I ) *DIRC0S (MEL- 1,1,2) 
0MEZ(INN+1)— G2WAKE(KCIR, J+l , I)*DIRC0S (MEL- 1 , 1,3) 
IF(M0D(NC0L, 2)  .EQ.  0)  THEN 
0MEX( INN+1 )-0.0D0 
OMEY ( INN+1 )-0.0D0 
0MEZ( INN+1 )-0.0D0 
END  IF 
END  IF 

1 40  CONTINUE 
150  CONTINUE 


o  n 


C 

C  Append  vorticities  for  the  midpoint  nodes  of  squares. 

★*■*■**•*** 


DO  180  I-l.NWPTS 

NCOL-NRWAKE(KCIR , I ) - 1 
DO  170  J-1.NC0L/2 
c 

c  Determine  node  #  of  midpoint 
c 

INM-NODE+(I-l)*(3*NCOL+2)+NCOL+l+J 

c 

c  Determine  node  #'s  of  4  corners  around  it 
c 

IN1-NODE+ ( I - 1 ) * ( 3*NC0L+ 2 ) +J 
IN2-IN1+1 

IN3-N0DE+ ( I - 1 ) * ( 3*NCOL+2 ) +2*NC0L+1+J 
IN4-IN3+1 

c 

c  Average  the  x- ,  y- ,  and  z-components  of  omega  vectors 
c  at  4  corners, 
c 

OMEX ( INM) - ( OMEX ( INI ) +0MEX ( IN2 ) +0MEX ( IN3 ) +OMEX ( IN4 ) ) /4 . 0  DO 
OMEY ( INM)— (OMEY ( INI ) +0MEY ( IN2)+OMEY( IN3 )+OMEY( IN4) ) /4 . 0D0 
OMEZ ( INM) - (OMEZ ( INI ) +0MEZ ( IN2 ) +OMEZ ( IN3 ) +OMEZ ( IN4 ) ) /4 . 0D0 
170  CONTINUE 
180  CONTINUE 
C 

C  Mirroring  vorticities  on  opposite  1/2  span  of  the  wake. 

C  NOTE:  for  symmetric  wings  only! 

C 

DO  200  I-l.NWPTS 

NCOL-NRWAKE (KCIR , I ) - 1 
DO  190  J-1.NC0L/2 

IN1-NODE+ ( I - 1 ) * ( 3*NCOL+2 ) +J 
IN2-NODE+ ( I - 1 ) * ( 3*NCOL+2 ) +2*NC0L+1+J 
INM-N0DE+( I - 1 ) * ( 3*NCOL+2 ) +NC0L+1+J 
K-NC0L-J+1 

IN01-N0DE+ ( I - 1 ) * ( 3*NCOL+2 )  +1+K 
IN02-N0DE+ ( I - 1 ) * ( 3*NCOL+2 ) +2*NCOL+2+K 
IN0M-N0DE+(I -l)*(3*NCOL+2)+NCOL+l+K 
OMEX(INOl)  — OMEX(INl) 

OMEY ( IN01 ) -OMEY ( INI ) 

OMEZ  ( INOl )  —  OMEZ  ( INI ) 

OMEX ( INO  2 )  —  OMEX ( IN2 ) 

OMEY ( IN02 ) -OMEY ( IN2 ) 

OMEZ ( IN02 )  —  OMEZ ( IN2 ) 

OMEX ( INOM ) - - OMEX ( I NM ) 

OMEY ( INOM ) -OMEY ( INM ) 

OMEZ ( INOM ) - - OMEZ (INM) 

190  CONTINUE 
200  CONTINUE 
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c  If  NCOL  is  odd,  append  missing  center  square  vorticities. 
c 

IF(MOD(NCOL, 2)  .EQ.  1)  THEN 

INL1-N0DE+ ( I - 1 ) * ( 3*NC0L+2 ) +NC0L/2+1 
OMEX ( INL1+ 1 )  —  OMEX ( INL1 ) 

OMEY ( INL1+1 ) -OMEY ( INL1 ) 

OMEZ ( INL1+1 )— OMEZ ( INL1 ) 

INL2-NODE+ ( I - 1 ) * ( 3*NCOL+2 ) + ( 5*NCOL) /2+2 
OMEX  ( INL2+1 )— OMEX  ( INL2 ) 

OMEY ( INL2 + 1 ) -OMEY ( I NL2 ) 

OMEZ  ( INL2+1 )— OMEZ  ( INL2 ) 

IMID-NODE+ ( I - 1 ) * ( 3*NCOL+2 ) + ( 3*NC0L) /2+2 
OMEX(IMID)— 0 . 0D0 
OMEY(IMID)— 0 . 0D0 
OMEZ ( IMID) -0 . 0D0 

END  IF 

C*************************************** 

C  Append  wake  c  factors  to  C FACT  array 

Q’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'kirk'k'k'k'k'k'k'k'k'k'k'k'k'k'kirk’kirk'k’k 

C 

c  Fill  with  ones 
c 

DO  220  I-NELE+1 , NWELE 
DO  210  J-1,3 

CFACT ( I , J ) -1 . 0D0 
210  CONTINUE 

220  CONTINUE 
c 

c  Fill  C-factor  array  for  wake 

c 

DO  240  I-NELE+1 .NWELE 
DO  230  J-1,3 
c 

c  Determine  node  number 
c 

NN-N ( I , J ) 
c 

c  Find  the  z-component  in  element's  local  frame  of  the  node's 
c  global  normal  (n-sub-z). 
c 

ANLZ-ANORX(NN)*DIRCOS(I , 3 , 1)+AN0RY(NN)*DIRC0S(I , 3,2) 
+ANORZ(NN)*DIRCOS(I ,3,3) 
c 

c  Convert  global  vorticity  at  node  (CAP  OMEGA)  into  element's 
c  local  frame  (little  omega), 
c 

0MLX-0MEX ( NN ) +DIRC0S (1,1,1) +OMEY ( NN ) +DIRCOS (1,1,2) 
+0MEZ(NN)*DIRC0S (I ,1,3) 

0MLY-0MEX(NN)*DIRC0S (1,2, 1)+0MEY(NN)*DIRC0S (1,2,2) 
+OMEZ(NN)*DIRCOS (1,2,3) 

OMLZ-OMEX(NN)*DIRCOS( I , 3 , 1)+0MEY(NN)*DIRC0S ( I , 3,2) 
+0MEZ(NN)*DIRC0S(I ,3,3) 
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c 

c  The  c-factor  is  calculated  using  Eq. (2.4-20)  of  Mracek's 
c  dissertation, 
c 

W  SQ— 0MLX*0MLX+0MLY*0MLY+0MLZ*0MLZ 
IF(WSQ  .GT.  1.0D-10)  THEN 

CFACT ( 1 , J ) -DSQRT ( WSQ ) /DSQRT ( WSQ+ ( 0MLZ*0MLZ) / ( ANLZ*ANLZ ) ) 
END  IF 

230  CONTINUE 
240  CONTINUE 

RETURN 

END 
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Appendix  C :  Subroutine  VELWK 

Subroutine  VELWK  calculates  the  disturbance  velocity  induced  at  a 
point  in  the  flow  by  the  wake  which  has  been  modeled  as  rows  of 
continuous  vortex  sheets.  The  routine  calls  the  already  existing 
routines,  VELBND,  for  each  triangular  panel  in  the  wake,  and  VELVF  for 
the  variable  strength  edge  cores  around  each  wake  row. 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  C 
c  ***************  c 
C  *  VELWK  *  C 
c  ***************  c 
c  c 

C  THIS  ROUTINE  CALCULATES  THE  VELOCITY  INDUCED  AT  A  POINT  C 
C  BY  THE  WAKE  WHICH  HAS  BEEN  MODELED  AS  ROWS  OF  CONTINUOUS  C 
C  VORTEX  SHEETS.  THE  METHOD  TRANSFORMS  GLOBAL  VORTICITY  AT  C 
C  A  NODE  INTO  LOCAL  COORDINATES  AND  MULTIPLIES  IT  BY  THE  C 
C  APPROPRIATE  VEL  OUTPUT.  C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


SUBROUTINE  VELWK 
IMPLICIT  REAL*8(A-H,0-Z) 

PARAMETER  (ME-1000 , MN-1000 , MEQ-200 , MUN-200 , MC-2 5 , MCI-100 ) 
PARAMETER  (MW-25 .MWP-25) 

COMMON/ELDATA/  NELE , NODE , N (ME , 3 ) ,XNODE(MN) , YNODE(MN) , ZNODE(MN) 
COMMON/EDGEDA/  NOPEN, NCLOSE,NTOTCR,NEDG(MC, MCI , 2) ,NCIRC(MC) 
COMMON/CONVEC/  NCOVCR,NCONCR(MCI) ,NCON(MC ,MCI , 2) .NCONPO(MCI) 
COMMON/CORNER/  AELE(ME) , BELE(ME) , CELE(ME) , DIRCOS (ME, 3 , 3) 
COMMON/VELINP/  XP , YP , ZP , IDENTE , NODE1 , NODE2 
COMMON /VELVFO/  VXVF(3) , VYVF(3) ,VZVF(3) 

COMMON/GNUMBE/  G1(MC,MCI) ,G2(MC,MCI) ,G3(MC,MCI) 

COMMON/GBYEDG/  GIW(MC.MCI) ,G2W(MC,MCI) ,G3W(MC,MCI) 

COMMON/ERRORS/  I ERROR , KERROR , MERROR , NERPOR 
COMMON/VELOUT/  VX(6) ,VY(6) ,VZ(6) 

COMMON/OMEGAL/  OMLOCX ( 3 ) , OMLOCY ( 3 ) 

COMMON/GNORMS/  ANORX(MN) , ANORY (MN) , ANORZ(MN) ,CFACT(ME,3) 
COMMON/CONSTS/  FOURPI , PI , CUTOFF , CUTWK 
COMMON/GOME GB/  OMEX(MN) , OMEY(MN) ,OMEZ(MN) 

COMMON/VCFINP/  XI , Y1 , Z1 , X2 , Y2 , Z2 , VXCF , VYCF , VZCF 
COMMON/STRWAK/  GAVE(MC.MCI) ,GWAKE(MC,MCI ,MW) , NRWAKE(MC , MW) 
COMMON/S PLITW/  ISPLIT(MC ,MCI ,MW) , DIMINI , ISPER(MCI ,MW) 
COMMON/TIMESD/  TIME , DTIME , NMOTIO , NCURTM , NTIME , NSTIME , NWPTS , MAXWK 
COMMON /VBOUND/  VBX.VBY.VBZ 
COMMON/VEWAKE/  VWKX , VWKY , VWKZ 
COMMON/CNTERS/  KCIR , NWELE 
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INITIALIZING  THE  VELOCITIES 

VWKX-0 . 0D0 
VWKY-0 . 0D0 
VWKZ-0 . 0D0 

LOOPING  THROUGH  THE  WAKE  ELEMENTS  FOR  EACH  ONES  CONTRIBUTION 

DO  10  I-NELE+1 , NWELE 
IDENTE-I 

FINDING  THE  LOCAL  VORTICITY  AT  EACH  NODE  OF  ELEMENT  I 
DO  20  J-1,3 

ANX-ANORX(N(I , J) )*DIRCOS (1,1,1) 
+ANORY(N(I,J))*DIRCOS(I,l,2) 
+AN0RZ(N(I,J))*DIRC0S(I,1,3) 

ANY-ANORX(N(I , J) )*DIRC0S(I ,2,1) 

+ANORY(N(I , J) )*DIRCOS(I ,2,2) 

+ANORZ(N(I , J) )*DIRC0S (I ,2,3) 

ANZ-ANORX(N(I , J) )*DIRCOS(I ,3,1) 

+ANORY(N(I , J) )*DIRCOS(I ,3,2) 

+ANORZ(N(I , J) )*DIRCOS(I ,3,3) 

Al 1-CFACT (I , J)*( ANX*ANX+ANZ*ANZ ) / ( ANZ*ANZ ) 
A12-CFACT ( I , J ) *ANX*ANY/ ( ANZ*ANZ) 

A2  2-CFACT (I , J) *(ANY*ANY+ANZ*ANZ ) / (ANZ*ANZ ) 
WOMEGX-OMEX (N ( I , J) ) *D IRCOS (1,1,1) 

+OMEY(N(I , J) )*DIRCOS(I ,1,2) 
+OMEZ(N(I,J))*DIRCOS(I,l,3) 

WOMEGY-OMEX (N ( I , J) ) *DIRCOS (1,2,1) 

+OMEY(N(I , J) )*DIRCOS (I ,2,2) 

+OMEZ(N( I , J) )*DIRCOS(I ,2,3) 

OMLOCX( J )-All*WOMEGX+Al2*WOMEGY 
OMLOCY ( J ) -Al 2*WOMEGX+A2  2  *WOMEGY 
20  CONTINUE 

CALLING  THE  VELOCITY  FOR  ELEMENT  I 
CALL  VELE 

ADDING  TO  THE  OTHER  ELEMENTS 
DO  30  J-1,3 

VWKX-VWKX+VX ( J ) *OMLOCX ( J ) 

+VX(J+3)*0ML0CY(J) 

VWKY-VWKY+VY ( J ) *0ML0CX ( J ) 

+VY(J+3)*0ML0CY(J) 

VWKZ-VWKZ+VZ ( J ) *0ML0CX ( J ) 

+VZ( J+3)*0ML0CY( J ) 

30  CONTINUE 

10  CONTIN^IE 
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c 

c 

c 

c 


c 

c 

c 

c 


c 


ADDING  THE  INFLUENCE  OF  THE  CORES  ALONG  THE  EDGES  OF 
EACH  WAKE  ROW. 

DO  50  I-NTOTCR+1 , NTOTCR+NWPTS 
DO  40  J-l.NCIRC(I) 

IF(J  .GT.  NCIRC(I)/2-I  .AND.  J  .LT.  NCIRC(I)- 
NODE1—NEDG ( I , J , 2 ) 

NODE2— NEDG ( I , J , 1 ) 

ELSE 

NODE I-NEDG ( I , J , 1 ) 

NODE2-NEDG ( I , J , 2 ) 

END  IF 

CALL  VELVF 

VWKX-VWKX+G 1 W ( I , J ) *VXVF ( I ) 

+G2W(I , J)*VXVF(2) 

+G3W( I , J ) *VXVF ( 3 ) 
VWKY-VWKY+G1W(I , J)*VYVF(1) 

+G2W(I , J)*VYVF(2) 

+G3W(I , J)*VYVF(3) 

VWKZ— VWKZ+G 1W ( I , J ) *VZVF ( 1 ) 
+G2W(I,J)*VZVF(2) 

+G3W( I , J ) *VZVF ( 3 ) 

40  CONTINUE 

50  CONTINUE 


RETURN 

END 


THEN 
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